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Executive  Summary 

The  goal  of  the  research  was  to  develop  a  new  methodology  for  the  fabrication  and  the  design  of  new 
multifunctional  composites  and  devices  using  3D  printing.  The  main  accomplishments  of  this  project  can 
be  summarized  as  follows:  1)  Proposing  and  demonstrating  the  4D  printing  concept  and  novel  printed 
active  composites  (PACs);  2)  developing  multiple  PAC  continuum  mechanics  PAC  models  of  various  levels 
of  accuracy  and  complexity;  3)  developing  the  4D  printing  concept  to  design  PACs  undergoing  multiple 
shape  changes;  4)  developing  a  PAC-based  concept  that  enables  reversible  shape  changes;  5)  application 
of  PACs  for  printed  origami  design;  6)  developing  a  finite  element  framework  for  the  optimum  design  of 
PACS  by  topology  optimization;  7)  optimizing  and  experimentally  characterizing  PACs  that  undergo 
desired  shape  changes. 

The  outcomes  of  this  research  project  were  published  in  14  peer-reviewed  journal  papers.  Our  paper  on 
topology  optimization  for  PACs  was  selected  featured  article  in  the  ASME  Journal  of  Mechanical  Design. 
Our  work  on  4D  printing  and  PAM  was  reported  by  a  variety  of  media,  including  NPR  and  ABC,  Denver 
Post,  and  others. 

1.  Summary  of  Major  Research  Activities 

1.1  Printed  Active  Composites  and  4D  Printing 

Recent  advances  in  3D  printing  allow  the  precise  placement  of  multiple  materials  at  micrometer  resolution 
with  essentially  no  restrictions  on  the  geometric  complexity  of  the  spatial  arrangement.  Complex  3D  solids 
can  be  created  with  highly  non-regular  material  distributions  in  an  optimal  fashion.  In  this  project,  we 
exploit  these  advances  and  introduce  the  paradigm  of  printed  active  composites  (PACs).  We  directly  print 
a  composite  from  a  CAD  file  that  specifies  the  inhomogeneous  material  layout  at  the  lamina  and  laminate 
level.  This  process  has  considerable  design  freedom  to  enable  creation  of  composites  with  complex  and 
controllable  anisotropic  thermomechanical  behavior  via  prescribing  the  material  layout  of  active  and 
inactive  phases.  Our  PACs  are  soft  composites  consisting  of  glassy  polymer  inclusions,  such  as  fibers, 
reinforcing  an  elastomeric  matrix.  The  glassy  polymer  phase  exhibits  a  shape  memory  effect  [1,  2]  and 
we  use  this  to  create  active  soft  composites  where  the  glassy  polymer  phase  serves  as  a  switch  to  affect 
the  shape  memory  behavior  of  the  composite.  In  this  manner  we  imbue  the  active  composites  with 
intelligence  via  the  lamina  and  laminate  architecture  and  the  thermomechanical  training  process.  The 
active  motion  of  the  PAC  is  essential  to  the  4D  printing  concept,  which  is  in  the  spirit  of  the  recent 
developments  of  Tibbits  [3],  although  our  work  differs  significantly  in  terms  of  the  physical  phenomena  at 
play  and  the  emphasis  on  understanding  them.  Upon  the  initial  development  of  4D  printing  and  PAC,  we 
further  pushed  the  design  space  of  4D  printing  by  developing  three  new  types  of  designs:  1)  sequential 
folding  of  components  using  digital  shape  memory  polymers;  2)  multiple  shape  change  printed  active 
composites  where  multiple  types  of  shape  memory  fibers  are  used  to  achieve  activation  of  the  composites 
at  different  time  and  different  conditions;  3)  reversible  shape  changing  composites  where  we  use  both 
shape  memory  polymers  and  hydrogels  to  achieve  autonomous  shape  change  without  the  assistance  of 
external  forces. 

Our  major  efforts  on  modeling  and  experimental  verification  included:  1)  In  order  to  gain  in-depth 
understanding  of  interplay  between  different  design  parameters  of  PAM,  we  developed  a  suite  of 
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theoretical  models  that  vary  in  phenomena  captured,  accuracy,  and  complexity.  These  models  were 
compared  against  experiments.  2)  We  investigated  the  effects  of  oxygen  on  the  interface  properties  of 
the  3D  printed  materials;  3)  Guided  by  the  theoretical  analysis,  we  applied  the  4D  printing  concept  to 
create  PAC  that  deform  into  complicated  geometries,  such  as  4D  printed  origami.  4)  We  explored  how  to 
use  the  frontal  polymerization  method  to  create  folding  structures.  5)  We  studied  methods  of  integrating 
conductive  ink  with  3D/4D  printing  techniques  to  enhance  the  multi-functionality  of  PACs. 

1.2  Design  and  Optimization  Methodology 

One  of  the  primary  goals  of  this  project  was  the  development  and  application  of  a  design  optimization 
methodology  for  determining  the  precise  and  optimal  layout  of  functional  composites  realized  by  3D/4D 
printing.  State-of-the-art  multi-material  printers  allow  fabricating  highly  non-regular  (but  carefully 
prescribed)  spatial  arrangements  of  multiple  polymers.  To  take  advantage  of  this  capability,  we  developed 
two  novel  3D  topology  optimization  approaches  that  determine  the  precise  geometry  of  each  material 
phase  within  the  composite.  The  optimization  method  accounts  for  the  complex  physical  phenomena  of 
PACs  via  proper  mechanical  models  and  high-fidelity  finite  element  predictions. 

We  studied  two  optimization  approaches;  one  describes  the  geometry  of  the  material 
arrangement  via  the  spatial  distribution  of  material  parameters  at  the  macro  scale.  The  second  approach 
resolves  the  material  architecture  via  level  set  functions  that  describe  the  microstructural  layout.  The  first 
approach,  often  called  density  approach,  has  been  used  successfully  in  the  past  for  a  broad  range  of 
structural  and  material  design  problems  [4]  but  not  yet  for  the  design  of  3D  PACs.  We  have  studied  the 
density  approach  for  the  design  of  PACs  approximating  the  shape  memory  behavior  of  the  active  phase 
by  a  simple  eigenstrain  model.  The  optimized  material  layouts  were  then  printed  on  an  Object  260  printer 
(Stratasys,  Edina,  MN,  USA).  This  study  served  mainly  as  proof-of-concept  and  revealed  several  issues  of 
the  density  method. 

Albeit  robust  and  easy  to  implement  for  simple  kinematic  and  constitutive  models,  the  main 
disadvantages  of  the  density  approach  are  that  (i)  it  cannot  resolve  reliably  complex  nonlinear  mechanical 
phenomena  and  (ii)  it  requires  post-processing  to  extract  the  final  material  layout  which  then  can  be 
printed.  To  avoid  these  issues,  we  developed  and  integrated  a  Level  Set  Method  (LSM)  used  to  describe 
the  microstructural  geometry  and  a  generalized  version  of  the  extended  Finite  Element  Method  (XFEM) 
for  discretizing  the  mechanical  models  and  predicting  the  PAC  response. 

We  developed  a  simple,  computationally  efficient,  but  still  sufficiently  accurate  model  to  predict 
the  PAC  response  at  the  training  stage  and  after  activation.  We  compared  this  new  model  against  an  even 
simpler  model  that  relies  on  crude  approximations  of  the  training  strains  and  verified  our  new  model 
against  experimental  data.  We  integrated  the  simplified  PAC  models  into  our  LSM-XFEM  optimization 
framework  and  applied  it  to  optimizing  PACs  that  assume  a  target  shape  upon  activation.  The  optimized 
designs  were  printed  and  their  behavior  studied  experimentally.  The  optimization  study  showed  that  the 
proposed  analysis  and  optimization  methods  are  well  suited  to  design  PACs  that  exhibit  a  desired 
response. 

We  extended  our  PAC  model  to  account  for  finite  strains  at  the  training  stage  and  after  activation. 
This  work  required  extensive  research  on  the  stability  of  the  XFEM  formulation  of  finite  strains  models 
and  we  developed  novel  stabilization  approach.  We  integrated  the  finite  strain  PAC  model  into  our  LSM- 
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XFEM  optimization  framework  and  performed  exploratory  studies  on  optimizing  PACs  that  assume  a 

target  shape  upon  activation.  Furthermore,  we  studied  the  influence  of  imperfect  bonding  between 

material  phases  on  the  optimum  design  by  incorporating  a  cohesive  zone  model  into  the  mechanical 

model. 

1.3  Publications 

The  outcomes  and  findings  of  this  project  were  published  in  the  following  peer-reviewed  journal  papers. 

[1]  *Ge,  0..,  Qi,  H .J.,  Dunn,  M.L.,  2013,  Active  Materials  by  4D  Printing.  Applied  Physics  Letter,  103, 
131901. 

[2]  Ge,  Q.,  Dunn,  C.,  Qi,  H  J.,  Dunn,  M.L.,  2014.  Active  Origami  by  4D  Printing.  Smart  Materials  and 
Structures,  23, 094007-15. 

[3]  Yu,  K.,  Ge,  Q.,  Qi,  H.J.,  2014,  Reduced  Time  as  a  Unified  Parameter  Determining  Fixity  and  Free 
Recovery  of  Shape  Memory  Polymers.  Nature  Communication,  5:3066. 

[4]  Yu,  K.,  Ritchie,  A.,  Mao,  Y.,  Dunn,  M.L.,  Qi,  H.J.,  2015.  Controlled  Sequential  Shape  Changing 
Components  by  3D  Printing  of  Shape  Memory  Polymer  Multimaterials.  Procedia  IUTAM,  12,  193- 
203. 

[5]  **K.  Maute,  A.  Tkachuk,  J.  Wu,  FI.  J.  Qi,  Z.  Ding,  and  M.L.  Dunn.  2015.  Level  Set  Topology 
Optimization  of  Printed  Active  Composites.  ASME  Journal  of  Mechanical  Design,  137:  111704-1, 
2015. 

[6]  Mao,  Y.,  Yu,  K.,  Isakov,  M.,  Wu,  J.,  Dunn,  M.L.,  and  Qi,  H.J.,  2015.  Sequential  Self-Folding  Structures 
by  3D  Printed  Digital  Shape  Memory  Polymers,  Scientific  Reports,  5:13616. 

[7]  Wu,  J.,  Yuan,  C.,  Ding,  Z.,  Isakov,  M.,  Mao,  Y.,  Wang,  T.,  Dunn,  M.L.,  Qi,  H.J.  2016.  Multi-shape 
active  composites  by  3D  printing  of  digital  shape  memory  polymers.  Scientific  Reports  6:24224. 

[8]  Mao,  Y.,  Ding,  Z.,  Yuan,  C.,  Ai,  S.,  Isakov,  M.,  Wu,  J.,  Wang,  T.,  Dunn,  M.L.,  and  Qi,  FH.J.,  2016.  3D 
Printed  Reversible  Shape  Changing  Component  with  Stimuli  Responsive  Materials.  Scientific  Reports 
6:24761. 

[9]  Zhao,  Z,  Mu,  X.,  Wu,  J.,  Qi,  FH.J.,  Fang,  D.,  2016.  Effects  of  oxygen  on  interfacial  strength  of 
incremental  forming  of  materials  by  photopolymerization.  Extreme  Mechanics  Letters  9(1):108— 
118. 

[10]  Zhao,  Z,  Wu,  J.,  Mu,  X.,  Qi,  FH.J.,  Fang,  D.,  2016.  Origami  Folding  by  Frontal  Photopolymerization, 
submitted. 

[11]  Mu,  Q,  Dunn,  C.K.  Wang,  L.,  Dunn,  M.L.,  Qi,  FH.J.,  Wang,  T.,  2016.  Thermal  cure  effects  on 
electromechanical  properties  of  conductive  wires  by  direct  ink  write  for  4D  printing  and  soft 
machines,  submitted. 

[12]  Sharma,  A.,  Villanueva,  FI.,  Maute,  2016.  On  Shape  Sensitivities  with  Fleaviside-Enriched  XFEM, 
submitted  to  Structural  and  Multidisciplinary  Optimization. 

[13]  Behrou,  R.,  Lawry,  M.,  Maute,  K.,  2016.  Level  set  topology  optimization  of  structural  problems  with 
interface  cohesion,  submitted  to  Int.  J.  Num.  Meth.  Engrg. 

[14]  Lawry,  M.,  Maute,  K.,  2016.  Level  set  shape  and  topology  optimization  of  finite  strain  bilateral 
contact  problems,  submitted  to  Int.  J.  Num.  Meth.  Engrg. 

*  This  work  was  widely  reported  by  the  media,  including  NPR-online,  ABC,  Denver  Post,  Engineering  and 
Technology  Magazine,  NASA  Tech  Briefs-Online,  etc.  At  the  end  of  November  2013,  it  was  estimated 
that  4.5M  people  was  reached  through  these  reports. 

**  Featured  paper  in  special  issue  of  ASME  Journal  of  Mechanical  design  on  additive  manufacturing  in 
December  2015. 
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2.  Research  Efforts  and  Accomplishments 

2.1  Fundamental  Models  and  Experimental  Studies 

2.1.1  Printed  Active  Composites  (PACs)  by  4D  Printing 

We  realize  the  PAC  concept  via  the  processing  approach  shown  in  Figure  2.1.1-la.  Briefly,  we  design  the 
complete  3D  architecture  of  the  fibers  and  matrix  in  a  CAD  file  to  yield  specific  active  behavior  based  on 
understanding  developed  from  theory  and  simulation,  and  use  a  3D  multimaterial  polymer  printer  (Objet 
Connex  260,  Stratasys,  Edina,  MN,  USA)  to  print  the  composite.  We  design  the  matrix  to  be  an  elastomer 
and  the  fibers  to  be  a  glassy  polymer  with  tailored  thermomechanical,  including  shape-memory,  behavior. 
This  design  is  printed  directly  and  results  in  a  film  that  contains  matrix  and  fiber  material  at  an  in-plane 
resolution  of  32-64  pm,  depending  on  the  printing  resolution  (Figure  2.1.1-lb).  The  complete  composite 
architecture  is  then  realized  by  printing  multiple  film  layers  to  create  an  individual  lamina  and  ultimately 
multiple  layers  to  create  the  3D  laminate.  We  measured  the  storage  modulus  Es  of  the  of  the  matrix,  fiber, 
and  the  composites  in  x-direction  with  a  dynamic  mechanical  analyzer  (DMA)  in  uniaxial  loading  at  0.1  Hz 
as  a  function  of  temperature  from  60  °C  to  15  °C  (Figure  2.1.1-lc).  The  matrix  modulus  decreases  linearly 
with  temperature,  consistent  with  rubber  entropic  elasticity;  the  slight  increase  at  15  °C  is  because  this  is 
close  to  the  glass  transition  region  and  continued  decreases  in  temperature  result  in  large  increases  in  Es. 


Figure  2. 1.1-1.  (a)  Schematic  illustrating  the  PAC  printing  process.  The  inkjet  heads  move  horizontally 
above  the  tray  depositing  multimaterial  droplets  of  polymer  ink  at  prescribed  positions,  wiping  them  into 
a  smooth  film,  and  then  UV  photopolymerizing  the  film.  After  one  film  layer  is  completed,  the  tray  moves 
down  to  print  the  next  layer,  (b)  Schematic  of  a  PAC  lamina;  fibers  are  oriented  at  an  angle  6  from  x- 
direction  (the  loading  direction),  (c)  The  storage  modulus  Es  in  x-direction  of  the  matrix  (M),  fiber  (F),  and 
the  composites  from  60  °C  to  15  °C  ,  along  with  the  modulus  obtained  from  quasi-static  uniaxial  tensile 
tests  at  60°C  and  15°C  (circles),  (d)  Theoretical  estimates  of  the  relaxed  modulus  Er  versus  fiber 
orientation,  including  experiments  at  vf  =  0.28  (circles)  at  15  °C. 
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The  composites  show  an  increasing  modulus  as  temperature  decreases,  with  a  large  variation  centered  at 
about  40  °C,  which  is  associated  with  the  glass  transition  of  the  fibers.  Considerable  anisotropy  exists, 
especially  as  the  fibers  transition  to  the  glassy  state  and  stiffen.  Also  shown  as  circles  in  Figure  2.1.1-lc 
are  results  for  the  modulus  obtained  from  quasi-static  uniaxial  tensile  tests  at  60  °C  and  15  °C;  they  are  in 
good  agreement  with  the  DMA  results.  Since  the  relaxation  contributes  importantly  to  the  shape  memory 
behavior  of  the  composites,  we  plot  the  relaxed  modulus  (after  30  minutes)  for  each  as  a  function  of  fiber 
orientation  in  Figure  2.1.1-ld.  Also  shown  are  theoretical  estimates  (solid  lines)  based  on  a  recent  theory 
of  fiber-reinforced  elastomers  (Guo  et  al.  [5,  6],  Lopez-Pamies  et  al.[7],  deBotton  et  al.[8]).  The  theory 
describes  the  experiments  well  and  both  show  significant  anisotropy  in  the  response. 

We  design  PACs  that  exhibit  the  shape  memory  effect  over  the  range  from  TL  =  15  °C  to  TH  =  60 
°C.  To  exercise  the  shape  memory  effect  (Figure  2.1.1-2a)  we  deform  samples  at  TH  =  60  °C  which  is  above 
7g  of  the  fibers  (~35  °C),  and  then  maintain  the  applied  stress  on  the  sample  while  cooling  it  to  TL  =  15  °C 
at  a  rate  of  2  °C/min.  At  TL  the  matrix  is  still  in  the  rubbery  state  but  the  fibers  have  transformed  to  the 
glassy  state.  During  cooling,  the  strain  decreases  slightly,  due  to  thermal  strain  of  the  constituents.  After 
holding  the  sample  at  TL  for  a  time  sufficient  to  equilibrate  the  strain  (~5  min),  we  released  the  stress  and 
observed  elastic  unloading,  the  magnitude  of  which  depends  on  the  fiber  orientation.  At  this  stage  the 
composite  is  in  its  fixed  configuration  and  we  characterized  the  degree  of  fixity  by 
Rf  ={su  ~sT)/{sm  ~£t)>  where  sm  and  su  are  the  axial  strains  before  and  after  the  unloading  step,  and 
sT  is  the  thermal  strain  at  Ti.  Finally,  we  recover  the  strain  by  heating  the  unconstrained  sample  at  2 
°C/min  to  60  °C.  Figure  2.1.1-2b  shows  the  applied  temperature  and  stress  cycle  and  Figure  2.1.1-2c  shows 
the  strain  response  for  the  composites  as  well  as  for  the  fiber  and  matrix,  indicating  a  significant 
anisotropy  in  the  shape  memory  behavior  of  the  composites  as  a  function  of  fiber  orientation;  this  is 
illustrated  specifically  for  the  fixity  in  Figure  2.1.1-2d.  The  fixity  as  a  function  of  6  follows  the  behavior 
of  the  relaxed  modulus  Er,  exhibiting  a  slight  minimum  at  0  ~  45-60°.  Predictions  of  a  theoretical  model, 
based  on  combining  models  of  hyperelastic  composites  with  those  describing  the  shape  memory  behavior 
of  glassy  polymers  are  in  good  agreement  with  measurements  for  the  shape  memory  cycle  and  fixity 
(Figure  2.1.1-2d).  We  note  that  the  polymer  fibers  themselves  exhibit  shape  memory  (  Rf  ~  80%),  while 

the  matrix  does  not  ( Rf  =  0%).  Simulations  show  that  fixity  increases  with  increasing  fiber  volume  fraction 

as  this  provides  increasing  stiffness  to  hold  the  composite  in  its  prescribed  deformation  state  after  cooling. 
For  all  volume  fractions,  the  fibers  are  most  effective  in  terms  of  providing  fixity  when  aligned  with  the 
load,  but  this  is  at  the  expense  of  reduced  strain  for  a  fixed  stress. 

We  further  exploit  the  anisotropic  shape  memory  behavior  of  the  active  lamina  to  create  PAC 
laminates  in  Figure  2. 1.1-3.  We  design  laminates  with  varying  fiber  orientations  and  volume  fractions, 
and  we  even  allow  these  design  parameters  to  vary  within  a  lamina,  an  aspect  that  can  be  readily  achieved 
via  printing.  Figure  2.1.1-3a  shows  the  notional  concept:  a  two-layer  laminate  with  a  prescribed  fiber 
architecture  is  printed,  then  heated,  stretched,  cooled,  and  released.  Upon  release  of  the  stress  it 
assumes  a  complex  temporary  shape,  depending  on  the  fiber  architecture.  Upon  reheating  it  recovers  its 
original  shape,  a  flat  rectangular  strip.  Figure  2.1.1-3b  shows  an  actual  strip  in  its  original  shape  and  then 
Figure  2.1.1-3c-h  shows  the  results  of  this  process  with  strips  that  are  nominally  identical  macroscopically, 
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At  7h,  both  the  matrix  and  fibers  are  soft. 
At  7i,  the  matrix  is  soft,  but  fibers  are  stiff. 


80  100 
Time  (min) 


0  15  30  45  60  75  90 

Fiber  Orientation  (°) 


Figure  2. 1.1-2.  Shape  memory  behavior  of  active  composite  lamina,  (a)  Schematic  of  a  shape  memory 
cycle,  (b)  Thermomechanical  loading  program  for  the  shape  memory  behavior  of  the  active  composite 
lamina,  (c)  Strain-time  response  of  the  vf  -  0.28  composites  (the  inset  denotes  the  matrix  (M)  and  fibers 
(F)).  (d)  Predictions  of  fixity  versus  fiber  orientation  are  compared  to  experiments  (circles). 


but  have  differing  fiber  architectures.  The  laminates  in  Figure  2.1.1-3c-e  each  have  two  layers  of  equal 
thickness  (0.5  mm),  one  of  them  has  a  fiber  volume  fraction  of  vf  =  0.25  and  the  other  is  pure  matrix  ( 

vf  =0),  and  6  -  0°,  90°,  and  30°  in  Figures  2.1. l-3c,  d,  and  e,  respectively.  Both  the  6  -  0°  and  90°  strips 

roll  into  coils,  with  the  6  -  0°  having  a  much  larger  curvature  due  to  the  increased  fixity  of  a  6  -  0°  lamina 
versus  a  6  -  90°  lamina.  In  general,  the  magnitude  of  the  curvature  depends  on  the  composite 
geometries  and  properties,  the  applied  mechanical  load,  and  the  thermal  history;  all  of  these  are  design 
variables.  The  strip  with  6  -  30°  coils  into  a  helix  due  to  the  shear-tension  coupling  of  the  laminate.  The 
laminates  of  Figure  2.1.1-3f,  g,  h  have  alternating  regions  of  segments  with  two  layers,  one  with  fibers  ( 
vf  -  0.25)  and  the  other  without  fibers  ( vf  =  0). 

The  ordering  of  the  segments  is  such  that  the  fibers  alternate  from  being  in  the  top  layer  and 
bottom  layer  along  the  strip.  In  Figure  2.1.1-3f,  the  fibers  are  at  0  =0°  and  so  the  strip  assumes  the 
shape  of  alternating  bent  segments,  while  in  Figure  2.1.1-3g,  the  fibers  are  at  6  =0°  in  the  top  layer  and 
at  6  -  30°  in  the  bottom  layer  leading  to  a  combination  of  bending  and  twisting  of  the  strip  in  alternating 
fashion.  In  Figure  2.1.1-3h,  the  fibers  are  all  at  6  =0°  but  are  patterned  in  segments  along  the  strip  of 
different  lengths  that  lead  to  a  variation  of  curvature  in  each  segment. 
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Figure  2. 1.1-4  shows  an  extension  of  this  concept  to  realize  a  shape  with  complex,  nonuniform 
curvature,  essentially  a  sculpted  surface.  Here  a  two-layer  laminate  (80  mm  x  80  mm)  is  fabricated  but 


Figure  2. 1.1-3.  Complex  low-temperature  shapes  of  active  composite  laminates  obtained  by  design  of  the 
laminate  architecture,  (a)  A  two-layer  laminate  designed  with  one  layer  being  a  lamina  with  fibers  at  a 
prescribed  orientation  and  one  layer  being  pure  matrix  material  is  printed,  then  heated,  stretched,  cooled, 
and  released.  Upon  release  of  the  stress  it  assumes  a  complex  shape,  depending  on  the  laminate 
architecture.  When  reheating  it  then  assumes  its  original  shape,  a  flat  rectangular  strip,  (b)  shows  an 
actual  strip  in  its  original  shape  and  (c)  -  (h)  show  results  of  this  process  with  differing  fiber  architectures. 


the  architecture  of  the  fibers  includes  a  spatial  variation  in  both  its  orientation  and  volume  fraction,  to 
yield  the  desired  curvature.  Over  each  end  of  the  laminate  0  =  0°  fibers  exist  at  a  volume  fraction  that 
varies  with  position  ranging  from  vf  =  0.25  at  the  end  of  the  laminate  to  vf  =  0.14  toward  the  center  of 

the  laminate  (the  precise  layout  is  shown  in  Supplementary  Material).  The  laminate  is  antisymmetric, 
with  respect  to  the  planar  dimensions.  The  surface  of  the  laminate  has  a  randomly  printed  patchwork, 
intended  to  be  reminiscent  of  the  interesting  artistic  and  architectural  forms  created  by  3D  printing  by 
Oxman  [9].  Unlike  that  complex  surface  that  was  directly  printed  in  3D,  ours  was  created  as  a  flat  plate 
and  then  following  our  thermomechanical  protocol,  it  assumed  the  curved  shape  as  a  result  of  the  fibrous 
architecture.  While  the  patchwork  pattern  here  is  only  for  aesthetic  purposes,  applications  can  be 
envisioned  where  the  pattern  has  functional  purposes. 
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Figure  2. 1.1-4.  A  sculpted  surface  with  a  complex,  nonuniform  curvature  obtained  by  design  of  the 
laminate  architecture,  (a)  Schematic  of  a  flat  laminate  that  is  stretched  at  TH.  (b)  After  cooling  to  71  and 
unloading,  a  desired  complex,  nonuniform  curvature  is  achieved. 


2.1.2  Continuum  Mechanics  PAC  Model 


Recently,  Guo  et  al.  [5,  6]  developed  a  hyperelastic  constitutive  model  for  anisotropic  fiber-reinforced 
composites,  which  treats  both  the  matrix  and  fibers  as  neo-Hookean  materials,  and  decomposes 
deformation  into  uniaxial  deformation  along  the  fiber  and  a  subsequent  shear  deformation  to  estimate 
the  strain  energy.  In  the  model,  the  strain  energy  stored  in  the  composite  system  is: 


fr  =  i//a(/1-3)  +  i/i/)(/4+/4-1/2-3),  (2. 1.2-1) 

and  the  corresponding  Cauchy  stress  is 

o  =  -pI  +  ^B  +  nb(l- 74“3/2) Fa0  ® Fa0 .  (2.1.2-2) 

iso 

aniso 


with  jua  = 


(1-v/)#«  +(1  +  v/)^/ 


and  nb=- 


v„vf 


-,  where  and  are  shear 


'(>  +  *,)«.  ”^(1  +  v/K+(1-v/)A/ . .  ~  " 

moduli  of  the  matrix  and  the  fibers,  respectively;  vm  and  vf  are  volume  fractions  of  the  matrix  and  the 


fibers  and  vm+vf  = 1.  Here,  F  is  the  deformation  gradient  F— Sx/cX  where  we  use  X  to  represent  the 

position  of  a  material  particle  in  the  original  (undeformed)  configuration,  while  x  is  the  position  of  the 
corresponding  particle  in  the  current  (deformed)  configuration.  B  is  the  left  Cauchy-Green  deformation 

gradient  and  B=FFT  .  C  is  the  right  green  Cauchy-Green  deformation  tensor  and  C  =  FTF.  Ix  and  I4 

are  the  principal  invariants  of  C:  7j=trC  and  IA  =  a0 -C-a0  =Ap,  where  a0  is  a  unit  vector  that 

represents  the  fiber  orientation  of  the  reinforced  fiber  in  the  original  (undeformed)  configuration  and  Ap 
is  the  fiber  stretch.  In  Eq.  (2),  p  is  an  arbitrary  hydrostatic  pressure,  the  second  term  is  associated  with 
the  isotropic  part  of  the  Cauchy  stress  and  the  third  term  is  associated  with  the  anisotropic  part  of  the 
Cauchy  stress.  Figure  lc  &d  show  the  comparison  between  model  simulations  and  experiments. 


Based  on  Guo's  model,  we  develop  a  thermomechanical  constitutive  model  to  describe  the  shape 
memory  behavior  of  composites  consisting  of  an  elastomer  matrix  and  glassy  polymer  fibers  that  pass 
through  the  glass  transition  during  the  temperature  cycles.  During  a  shape  memory  cycle,  the  model 
decomposes  the  total  deformation  F  into  the  thermal  deformation  Fr  and  the  mechanical  deformation 
Fm,  that  is,  F  =  FrFM  .  The  thermal  deformation,  Fr,  is  directly  described  by  a  simple  model  based  on 
the  observation  of  experiments.  The  mechanical  deformation,  FM ,  gives  rise  to  the  stress  acting  on  the 
composite  system.  At  high  temperatures,  fibers  are  in  a  rubbery  state,  and  //y  is  (this  value  is 
measured  by  experiments).  Therefore,  at  time  t  =  t0  with  a  high  temperature  T  =  T0 ,  the  Cauchy  stress  is 


<*.  -  -pi + P.U.,.,  f;  K  )t + Mb\  )F”: 


*  F,°,an 


(2. 1.2-3) 


During  glass  transition,  for  the  sake  of  model  simplicity  we  make  several  assumptions  and: 
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i) 

N) 


We  phenomenologically  assume  that  the  fibers  can  be  described  as  an  aggregate  of  rubbery  and  glassy 
phases,  and  the  phase  transition  between  these  two  phases  is  realized  through  the  change  of  volume 
fraction  of  each  phase  [10]. 


The  volume  fraction  of  each  phase  (  /'  for  the  rubbery  phase  and  f  for  the  glassy  phase)  as  a 
function  of  temperature  is  defined  as: 


fr  = 


i+exp[-(r-7;)/^] 


,  and  /„  =1  -fr, 


(2. 1.2-4) 


with  a  parameter  A  characterizing  the  width  of  the  phase  transition  zone,  and  a  reference 
temperature  Tr . 


iii)  During  glass  transition,  //,  is  simply  described  as  juf  -  fr/J.r+  fgjug  ( jug  is  also  measured  by 


experiments).  As  a  result nf  varies  from  pr  to  /Jg,  which  leads  to  the  change  of  jua  and  jub  in  Eq. 
(2),  from  jua  I  [jub I  )  to  //a  I  (jub I  ). 

~a  'Mf-Mr  ^b'Vf=Mr  a'i“f=Mg  ^b'Mf=Ms 

iv)  We  assume  that  the  newly  formed  glassy  phase  is  in  a  stress-free  configuration  and  take  this  as  its 
reference  configuration  for  any  subsequent  deformation  7 . 


Based  on  the  assumption  that  the  newly  formed  glassy  phase  is  in  a  stress-free  configuration,  when 
we  calculate  the  total  Cauchy  stress  applied  to  the  composite  during  glass  transition,  the  increasing 
amounts  of  pa  and  pb  is  only  associated  with  the  deformation  after  the  formation  of  the  glassy  phase 
resulting  in  the  corresponding  increase  of  /ua  and  /ub.  Therefore,  we  obtain  a  general  expression  for  the 
total  Cauchy  stress  acting  on  the  composite  during  glass  transition,  at  time  t-t0  +  nAt ,  and  T  -T0  +  nKT 


a„  =  -pi  +  (K  )T  +  Mbfi  [l  -  h  (F£,a0  f" 

+Z  { A^F T  (Fr  )T  +  A/4.,  [l  -  /4  (Fr ,  a,,  f" 


F>0®F>0 


F-X  0  F;^,. 


f  n 


In 


Eq.  (2. 1.2-5),  Ha(b)  Ma(b)\^{t=toy  AK(by  Ma&) |A/(t=<i+/At)  L/((=(o+(i_1)Ai) '  11^* 


M 

V  1=1  J 


(2. 1.2-5) 


F'O 

M  > 


riAFi 

V  J=i  J 


 F^ap 


,  a,  =  |— r — — r ,  and,  /4(F,a)  =Fa-Fa.  F®,  is  the  mechanical  deformation  before  glass 

rM  ao 


transition,  and  AF^  is  the  incremental  deformation  at  time  t  =  t0+iAt .  Comparing  model  predictions 
and  experiments  showed  good  agreement;  see  Figure  2. 1.1-2. 


2.1.3  Origami  by4D  Printing 

We  further  extend  the  PAM  concept  to  demonstrate  how  we  can  print  flat-plate  structures  consisting  of 
PAC  hinges  directly  connected  to  rigid  plastic  components  of  arbitrary  shape  and  then  program  the  hinges 
to  assemble  the  as-printed  structure  into  a  desired  3-D  configuration.  First  we  designed  a  box  consisting 
of  six  sides  connected  by  PAC  hinges  that  is  printed  in  a  flat  (unfolded)  form  as  shown  in  Figure  2.1.3-la. 
Specifically,  we  created  hinges  from  PACs  by  following  a  theoretical  study  similar  to  the  one  in  Section 
2.1.2,  and  stretched  the  box  biaxially  by  20%.  Figure  2.1.3-la  shows  the  as-printed  box  (the  flat  plate) 
where  the  rigid  sides  are  white  and  the  hinges  are  black.  The  assembled  box  was  created  after  biaxially 
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stretching  the  as-printed  structure  by  20%  at  TH ,  cooling  to  TL,  and  releasing  the  load.  It  clearly 
assembles  into  the  desired  box  shape  (Figure  2.1.3-lb),  with  only  small  deviations  from  the  desired  90° 


(a) 

f*0 

(b) 

/  €° 

o 

A _ 

*  Heating  1  | 

{*0 

(C) 

■A- 

f  £o 

(d) 

*0 

O 

4  Heating 

l£# 

Figure  2. 1.3-1.  Active  origami  box  and  pyramid.  The  printed  flat  cross  shape  in  (a)  assembles  itself  into  a 
desired  box  shape  in  (b)  after  the  programming  steps.  The  printed  flat  Ninja  star  shape  plate  in  (c) 
assembles  itself  into  a  desired  pyramid  shape  in  (d)  after  the  programming  steps. 

angles,  and  these  are  likely  due  to  inaccuracies  in  the  straining  process.  Figures  2.1.3-lc  and  d  show  a 
similar  structure,  a  five-sided  3D  pyramid.  Flere  the  printed  flat  Ninja  star  shape  plate  consists  of  a  square 
base  with  four  triangular  sides  (Figure  2.1.3-lc)  that  are  folded  to  a  3D  pyramids  with  60°  angles  (Figure 
2.1.3-ld). 

We  can  also  create  complex  3D  shapes  with  different  hinge  angles  by  printing  hinges  with 
different  geometries.  Figure  2.1.3-2a  shows  an  example  that  a  flat  triangle  sheet  that  folds  itself  into  an 
origami  airplane  with  a  0°  angle  in  the  middle  hinge  that  bends  upward  and  90°  angles  in  the  two  side 
hinges  that  bend  downward.  The  hinges  are  created  by  printing  the  bilayer  composites  with  the  fibers  on 
the  top  layer  for  the  0°  hinges  and  on  the  bottom  for  the  90°  hinges.  To  simplify  the  loading  process 
(allowing  simply  a  20%  stretch)  we  print  the  90°  hinges  in  Figure  2.1.3-2a  at  an  inclined  angle  relative  to 
the  base  of  the  plane,  but  we  maintain  the  fiber  orientations  parallel  to  the  base  (and  the  applied  stretch, 
Figure  2.1.3-2a  inset).  We  create  an  even  more  sophisticated  origami  airplane  by  printing  not  only  hinges 
with  different  lengths  but  also  ones  with  different  cross-section  profiles  (Figure  2.1.3-2c  and  d). 

The  use  3D  printing  of  active  materials  to  create  components  that  controllably  change  their  shape  over 
time  is  not  limited  to  the  use  of  PAC  hinges.  In  fact,  we  can  directly  print  3D  devices  by  strategically  placing 
shape  memory  polymers  at  pivotal  locations  or  throughout  an  entire  structure.  We  can  then  program  a 
temporary  shape  of  arbitrary  form  that  can  be  achieved  by  applying  a  prescribed  mechanical  loading  at 
Th  followed  by  releasing  the  constraints  at  TL  .  The  components  can  then  be  returned  to  their  complex 
original  3D  shapes  after  heating  back  to  TH.  Figure  2. 1.3-3  shows  such  an  example.  We  directly  printed  a 
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3D  box  with  a  pattern  of  SUTD-CU  logos  on  the  five  panels.  In  Figure  2.1.3-3a,  the  complete  3D  box  was 
printed  with  a  three-layer  laminate  consisting  of  a  Verowhite  middle  layer  embedded  in  a  Tangoblack 


Figure  2. 1.3-2.  Active  origami  airplanes.  A  flat  triangle  sheet  with  three  hinges  in  (a)  assembles  itself  into 
an  origami  airplane  with  a  0°  angle  in  the  middle  hinge  that  bends  upward  and  90°  angles  in  the  two  side 
hinges  that  bend  downward  in  (b).  A  flat  triangle  sheet  with  five  hinges  in  (c)  assembles  itself  into  an 
origami  airplane  with  two  winglets  in  (d). 


matrix.  The  box  is  then  deformed  to  a  flat  form  (Figure  2.1.3-3b)  by  applying  mechanical  loads  at  TH,  and 
cooled  to  Tl  where  the  constraints  were  removed,  leaving  it  in  a  flat  form.  Upon  heating  back  to  TH,  the 
structure  retakes  the  original  3D  box  shape  (Figure  2.1.3-3a). 
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Figure  2. 1.3-3.  A  directly-printed  origami  SUTD-CU  box.  An  as-printed  3D  SUTD-CU  box  in  (a)  was 
deformed  into  a  flat  form  at  TL  in  (b).  After  heating  back  to  TH,  the  structure  recovers  the  3D  box  shape. 

2.2.1  Sequential  folding  structures  using  digital  shape  memory  polymers 

In  this  work,  we  used  digital  SMPs  with  different  thermomechanical  properties  to  achieve  different 
recovery  speeds  of  the  SMPs  at  the  same  recovery  temperature;  we  then  further  exploited  this  feature  to 
design  complicated  folding  structures.  This  work  was  published  in  "Mao,  Y.,  Yu,  K.,  Isakov,  MS.,  Wu,  J.,  Dunn , 
M.L,  Qi,  H.J. ,  2015.  Sequential  Self-Folding  Structures  by  3D  Printed  Digital  Shape  Memory  Polymers,  Scientific 
Reports,  5:13616.  doi:10.1038/srep!3616." 


Hinge  Folding:  Hinges  designed  with  digital  SMPs  with  specified  thermomechanical  behaviors,  especially 
recovery,  are  the  key  element  of  sequentially  self-folding  structures  (shown  in  the  inset  of  Figure  la).  We 
designed,  fabricated,  and  tested  hinges  using  seven  different  digital  SMPs.  We  use  H-i  (i=l,2,3...7)  to 
represent  the  /- th  hinge.  The  glass  transition  temperatures  (Tgs)  of  these  hinge  materials  are  shown  in 
Figure  la.  The  shape  memory  behavior  of  the  hinges,  including  shape  fixity  and  shape  recovery,  is 
controlled  by  the  applied  mechanical  and  thermal-temporal  programming  and  recovery  conditions.  To 
understand  the  basic  hinge  behavior  we  consider  a  hinge  printed  with  an  initial  angle  of  0  =  90°  and  then 
programmed  (by  bending  at  TH  =90  °C)  to  assume  a  temporary  flat  state  ( 0  =  180° ),  which  corresponds 
to  a  shape  fixity  of  100%.  Then  we  recover  the  hinges  to  their  original  shape  ( 0  =  90° )  by  heating  them 
in  a  constant  temperature  bath,  during  which  we  monitor  the  recovery  of  the  hinges  as  a  function  of  time. 
We  quantify  the  shape  recovery  with  the  angular  recovery  ratio  (Rr): 


m= 


180°  -0(7) 
90° 


(2.2. 1-1) 


where  6(t)  is  the  time  dependent  angle  (Figure  2.2.1-la). 

The  time-dependent  angular  recovery  of  the  bending  hinges  can  be  examined  by  using  a  simple 
scaling  approach  with  respect  to  the  uniaxial  shape  memory  behavior,  which  is  modeled  by  a  multi-branch 
constitutive  model  4345  with  16  branches.  Figure  lb  shows  the  folding  time  at  three  different  folding 
temperatures  for  the  seven  hinges  by  ID  model.  Flere,  the  folding  time  is  defined  as  the  time  needed  for 
a  hinge  to  fold  to  0  =  91.8° ,  which  corresponds  to  a  recovery  ratio  of  98%.  It  can  be  seen  that  at  the  same 
folding  temperature,  hinges  fabricated  by  lower  Tg  materials  (the  higher  hinge  numbers)  need  longer  time 
to  fold.  The  difference  in  the  folding  time  provides  the  opportunity  to  design  the  sequencing  of  folding. 
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Figure  2.2. 1-1.  (a)  Glass  transition  temperatures  of  the  seven  digital  SMPs  used  for  hinge  materials;  the 
inset  shows  a  schematic  view  of  the  folding  of  an  SMP  hinge;  (b)  folding  time  of  hinges  to  obtain  98% 
recovery  ratio  at  various  recovery  temperatures. 

Sequential  self-folding.  As  the  first  example,  we  demonstrate  a  seemingly-simple  helical  structure  where 
rigid  (non-active)  panels  are  connected  by  active  hinges  on  each  corner  with  a  radius  of  5mm  (Figure 
2.2.1-2a).  The  thickness  of  the  panels  and  hinges  is  0.8mm,  the  depth  is  6mm  and  the  clearance  between 
neighboring  panels  is  5mm.  We  design  two  structures  that  have  identical  geometry,  but  each  has  a 
different  set  of  materials  used  for  the  hinges.  In  Sample  1  each  of  the  9  hinges  uses  the  same  SMP  material 
(H-3),  and  in  Sample  2  the  hinges  are  made  from  seven  different  materials  (H1-H7;  Table  2. 2. 1-1).  After 
printing  the  structure,  we  deform  it  into  a  flat  configuration  in  hot  water  at  TH= 90  °C,  which  is  above  the 
Tg  of  all  of  the  SMP  sections.  Then  the  sample  is  cooled  to  TL=10  °C,  at  which  all  the  SMP  hinges  are  in 
their  glassy  states.  After  releasing  the  external  load,  the  structure  is  fixed  at  the  temporary  straight  shape. 
To  activate  the  shape  recovery  of  the  structure,  we  immerse  it  in  hot  water  with  Tr= 90  °C  so  that  all  hinges 
experience  the  same  thermal  conditions.  The  shape  recovery  process  is  monitored  by  a  video  camera. 

Figure  2.2.1-2b  shows  the  recovery  sequence  of  Sample  1.  Since  all  the  hinges  possess  the  same 
thermomechanical  behavior  (H3  with  Tg=  ~55  °C),  they  start  to  recover  their  shape  simultaneously.  The 
uncoordinated  and  simultaneous  motions  of  all  hinges  lead  to  two  undesirable  results.  First,  at  ~7.5s  parts 
of  the  structure  collide  with  each  other  because  of  folding  pathway  interference.  The  contact  between 
different  parts  of  the  structure  provides  frictional  resistance  and  leads  to  a  significant  reduction  of  folding 
speed.  Second,  at  12.5s,  the  frictional  force  becomes  so  large  that  folding  stops  and  the  structure  fails  to 
reach  its  final  target  shape. 

Table  2.2. 1-1  Hinge  assignment  of  each  fold 

Fold  Number  1  2  3  4  5  6  7  8  9 

.  Sample  1  H-3  H-3  H-3  H-3  H-3  H-3  H-3  H-3  H-3 

atena  Sample  2  H-l  H-l  H-2  H-3  H-4  H-5  H-6  H-7  H-7 


For  the  structure  shown  in  Figure  2.2.1-2a,  in  order  to  recover  its  final  target  shape,  the  inner  hinges 
require  faster  folding  speeds  than  those  of  the  outer  hinges,  which  can  be  accomplished  by  assigning  SMPs 
with  lower  Tg s.  To  this  end  we  designed  hinges  unique  to  each  corner  as  (see  Case  2  in  Table  2.2. 1-1). 
After  experiencing  the  same  programming  procedure  as  described  in  the  previous  section,  the  structure 
is  immersed  in  hot  water  at  TH=  90  °C  and  the  free  recovery  process  is  monitored.  Figure  2.2.1-2c  shows 
the  snapshots  at  different  recovery  times,  which  visually  demonstrates  the  sequential  recovery  of  hinges 
to  yield  unconstrained  and  rapid  (~7.0  s)  folding  to  the  target  configuration  as  shown  in  Figure  2.2.1-2c. 
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Figure  2.2. 1-2.  Schematic  of  the  helical  SMP  component  (a).  Series  of  photographs  showing  the  shape 
recovery  process  of  the  helical  SMP  component  with  uniform  hinge  section  (b),  and  with  graded  hinge 
sections(c). 

The  different  folding  speeds  of  the  hinge  sections  due  to  their  respective  Tg s  enable  the  avoidance  of  self¬ 
collisions.  Since  the  folding  speeds  of  the  three  inner  hinges  are  faster  than  the  outside  ones,  they  first 
exhibit  shape  recovery  (they  start  the  shape  recovery  within  ~ls  and  finish  in  ~3s).  Then  shape  change  is 
successively  triggered  along  the  helical  line  and  the  last  two  hinges  start  to  recover  after  3s  and  finish  at 
~7.0s.  The  hierarchical  shape  recovery  profiles  of  the  hinge  sections  enable  the  successful  sequential 
shape  recovery  in  the  SMP  component  without  interruption  by  collisions. 

The  above  two  cases  demonstrate  the  importance  of  design  and  the  ability  to  control  the 
assembly  sequence  and  dynamics  by  the  spatial  arrangement  of  SMPs  with  varying  thermomechanical 
behavior.  A  "bad"  design  without  careful  coordination  of  different  folding  panels  can  lead  to  slower  and 
even  failed  folding  of  the  entire  structure. 

Folding  simulations:  finite  element  analysis  and  reduced-order  model.  To  describe  both  unfolding  and 
folding  behaviors  and  use  this  understanding  in  current  design,  we  take  a  two-stage  approach  for  analysis. 
First,  we  simulate  the  folding  dynamics  using  finite  element  analysis  (FEA).  FEA  provides  high-fidelity 
simulation  of  the  spatial  and  temporal  variation  of  the  programming  (unfolding)  and  deployment  (folding) 
processes.  Flowever,  it  is  also  time-consuming  as  it  involves  complicated  contact  (used  for  unfolding)  as 
well  as  nonlinear  geometrical  mechanics.  Second,  we  develop  a  reduced-order  model  (ROM)  to  simulate 
rigid  body  motions  of  the  panels  (driven  by  the  recovery  of  angular  shape  recovery  of  hinges)  to 
investigate  the  folding  pathways.  During  folding  and  unfolding  the  deformation  is  mainly  localized  in  the 
hinges  so  in  the  ROM  we  idealize  the  structure  as  rigid  panels  connected  by  hinges  with  rotational 
behavior  given  by  our  constitutive  model. 

For  ROM,  we  use  the  scaling  rule  demonstrated  above  to  predict  the  shape  recovery  ratio  of  hinge 
by  the  recovery  behavior  of  a  sample  under  linear  stretch.  Flere,  we  assume  that  the  folding  structure  can 
be  represented  by  rigid  panels  connected  by  hinges,  which  occupy  infinitesimal  space  but  control  the 
angles  between  two  panels.  Figure  2.2. 1-3  shows  panels  are  connected  by  hinges,  which  are  represented 
by  Pj  with  its  coordinates  of  (Xj,yi)  for  the  i-th  hinge.  In  the  ROM,  we  assume  that  the  angular  rotation  of 
the  hinge  is  equal  to  the  linear  recovery  ratio  of  a  SMP.  This  assumption  is  validated  in  the  above  sections. 
Following  Eq.  2.2. 1-1,  we  have 

e.(t)=mo-%oRi(t),  (2.2.1-2) 

where  6.{t )  is  the  angle  between  the  i-th  panel  and  (/+l)-th  panel,  R!  (t )  is  the  recovery  ratio  of  the  /- th 
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Figure  2.2. 1-3.  Schematic  of  two  arbitrary  strips  having  a  factitious  intersection  point 

hinge,  which  is  between  the  /- th  and  (/+l)-th  panels.  In  order  to  calculate  the  location  of  the  hinges,  the 
orientation  of  the  i-th  panel  is  calculated  as 

=  +  (2.2.1-3) 

where  ^.(7)  is  the  orientation  of  /-th  panel.  Finally,  the  coordinates  of  the  /-th  hinges  are 

i  i 

xi(t)  =  xo  +Z/;C0SW)]'  ^,(0  =  +  (2.2.1-4) 

i=l  i=\ 

where  /  is  the  length  of  the  /-th  panel. 

In  ROM  simulations,  at  each  time  increment,  the  coordinates  of  individual  hinge  points  are  first 
calculated;  the  intersections  points  and  the  intersection  indices  are  calculated.  If  0  <  Sg  <  1 7  the  above 

three  conditions  will  be  checked;  if  any  one  of  them  is  satisfied,  there  is  no  intersection;  otherwise,  an 
intersection  exits  and  folding  cannot  be  completed. 

Figure  2.2.1-4a  shows  finite  element  simulations  of  the  folding  and  unfolding  processes.  In  the 
unfolding  process,  the  structure  is  opened  by  rigid  surfaces  to  successively  unfold  one  panel  after  the 
other  at  a  high  temperature.  After  lowering  the  temperature  and  unloading,  the  structure  is  fixed  into  the 
flat  configuration.  Recovery  is  achieved  by  heating  the  structure  uniformly  and  it  folds  into  the  original 
configuration.  The  simulated  recovery  ratios  for  representative  hinges  (H-l,  H-4,  and  H-7)  are  presented 
in  Figure  2.2.1-4b  along  with  measurements  and  results  from  the  ROM.  Both  FEA  and  ROM  predictions 
are  in  good  agreement  with  experimental  measurement.  Figures  2.2.1-4c  and  2.2.1-4d  show  ROM 
simulations  of  the  two  samples  in  Table  2.2. 1-1  .  They  describe  both  the  unsuccessful  (Fig.  2.2.1-4c)  and 
successful  (Fig.  2.2.1-4d)  folding  sequences  well.  Importantly  the  ROM  captures  the  self-collisions  that 
occur  during  the  failed  folding;  based  on  this  capability,  we  use  it  to  design  more  complex  folding  scenarios 
in  the  next  section. 

Collision  determination.  Utilizing  the  ROM  simulation,  we  develop  a  simple  collision  index,  which  can  be 
used  to  determine  if  two  panels  will  collide  during  the  folding  process.  We  consider  two  arbitrary  panels 
in  the  structure  (as  shown  in  Figure  2.2. 1-3),  denoted  as  with  coordinates  of  points  as 

Pf]  (xf ,  ) #  p(j2)  (xf] ,  yf] ) #  respectively.  The  parametric  equation  for  the  /-th  panel  (denoted  as  a 

solid  line)  is  given  as 
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where  S  is  a  parameter  that  is  the  distance  from  point  /  to  point  (x,y)  normalized  by  the  length  of  the 
panel.  Therefore,  for  the  point  (x,y)  on  the  panel,  0<S<1.  By  using  the  parametric  equations  of  the  /-th 
panel  and  the  y'-th  panel,  the  coordinates  of  the  intersection  between  these  two  panels  can  be  calculated 

as  Cy  ( Xjj^jj ) .  Substituting  the  coordinate  of  intersection  point  (C9  <-W)  into  the  /-th  line  equation 


(Eq.(  2.2. 1-2))  yields  a  parametric  index  St .  When  0  <  St  <  1 ,  the  intersection  point  is  on  the  /-th  panel, 
while  when  S[  >  1  or  Sf  <  0  the  intersection  point  sits  on  the  extension  line  of  the  /-th  panel.  The  same  is 
true  for  parametric  index,  Sj ,  for  the  y’-th  panel.  It  is  thus  easy  to  conclude  that  when  the  two  panels 
intersect,  we  have  0  <  St  <  1  and  0<  Sj<  \  ■  otherwise  there  is  no  intersection. 

Closed  Shape  r^=:r  90  " 


(c)  (d) 

Figure  2.2. 1-4.  (a)  Finite  element  simulation  of  shape  memory  recovery  activated  sequential  self-folding 
when  the  strip  is  firstly  programmed  at  high  temperature,  (b)  Simulation  of  folding  angles  of  three 
representative  hinges  by  ROM  and  compared  with  that  by  FEA  and  experiments.  Solid-line  only:  FEA; 
round  dot:  ROM;  square:  experiments.  Red  color:  H-l;  magenta  color:  FI-4;  Blue  color:  FI-7.  Comparison 
between  ROM  simulations  and  experiments  for  Sample  1(c)  and  2  (d),  where  red  lines  indicate  the 
simulation  results. 
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We  further  define  an  interaction  index  for  intersection  of  line  i  and  line  j  as 

Sij  =  Si  x  Sj  .  (2.2.1-6) 

In  general,  when  0  <  Sy  <  1 ,  intersection  occurs.  However,  three  cases  that  can  give  0<Sy  <  1  but  no 
intersection:  (1)  0  <  S;<  1  and  Sj  >1  (for  example,  5;  is  slightly  larger  than  1);  (2)  0  <  Sj  <  1  andS,>l(for 

example,  5,  is  slightly  larger  than  1);  (3)  Si<  0  and  Sj  <0.  As  described  in  the  Methods  section,  the 
collision  index,  S#y,  is  calculated  based  on  the  motion  of  any  pair  of  panels  (the  /-  th  panel  and  the  j-  th  panel) 
during  each  time  increment  of  the  simulation.  In  general,  when  0  <  Sy  <  1 ,  a  collision  will  occur.  Figures 

2.2.1- 5a  and  2.2.1-5b  show  the  collision  indices  of  the  two  designs  discussed  above;  there  are  10  panels 
in  each  design  which  gives  rise  to  40  unique  panel  pairs.  In  Figure  2.2.1-5a,  the  structure  can  fold  without 
collisions,  therefore,  no  collision  index  falls  between  the  band  between  0  and  1  during  the  simulation.  In 
Figure  2.2.1-2c  two  collisions  occurred  between  panels  and  the  structure  cannot  fold  to  its  target  shape; 
Figure  2.2.1-5b  reflects  this  as  two  collision  indices  fall  in  span(0,l).  In  the  figure,  the  red  dot  between 
(0,1)  corresponds  to  first  collision;  the  blue  dot  corresponds  to  the  second  collision.  More  collision  will 
happen;  but  since  our  interest  is  the  first  collision,  no  further  calculation  was  conducted. 

A  self-locking  system.  Here  we  use  our  ROM  to  design  and  then  3D  print  a  self-folding  and  self-locking 
structure,  as  shown  in  Figure  6a.  Here,  two  holes  with  different  dimensions  are  designed  on  one  of  the 
end  panels.  The  objective  is  to  sequentially  guide  the  other  end  of  the  structure  through  the  first  hole  and 
then  through  the  second  hole  to  lock  the  structure  in  place.  Five  hinges  (with  a  radius  of  5mm  in  uniform) 
are  used.  In  order  to  ensure  that  the  end  of  the  strip  passes  through  the  two  holes  in  sequence,  we  design 
each  hinge  to  control  the  folding  rate  of  each  segment  of  the  structure.  This  results  in  two  requirements: 
First,  the  folding  of  hinges  1  and  2  should  be  faster  than  the  others  so  that  the  two  end  sides  reach  each 
other  in  the  correct  position  and  angle.  Second,  folding  of  hinges  4  and  5  should  be  slower  to  guarantee 
them  to  pass  through  the  two  holes  successively  and  then  lock  the  structure.  Using  our  ROM  simulations 
with  the  collision  index,  we  determined  feasible  folding  rates  for  each  hinge,  and  then  realized  them  with 
hinges.  Specifically,  we  assigned  hinges  H-2,  H-2,  H-3,  H-6  and  H-6  for  hinges  1  to  5  in  Fig.  6a.  The  self¬ 
folding  and  self-locking  structure  is  demonstrated  experimentally  in  Fig.  5.  The  structure  shown  in  Figure 

2.2.1- 6a  is  first  printed.  At  the  temperature  of  90  °C,  it  is  opened  into  a  flat  strip,  or  namely  the  unlocked 
configuration.  It  is  then  cooled  to  10  °C  where  the  flat  shape  is  maintained.  After  this  programming  step, 
the  unlocked  structure  is  immersed  in  hot  water  (90  °C)  and  the  free  recovery  process  is  shown  in  the 
snapshots  of  Figure  2.2.1-6d.  The  insets  show  the  ROM  simulated  folding  and  locking  procedure.  As 
predicted  by  ROM,  the  active  locking  side  is  able  to  pass  through  the  two  holes  successively,  and  finally 
recovers  to  the  initial  configuration  as  shown  in  Figure  2.2.1-6c.  The  ROM  simulation  and  the  experiment 
show  good  agreement. 

A  3D  folding  structure.  Folding  boxes  are  notable  folding  structures  in  daily  life  and  widely  used  in 
packaging.  They  normally  require  a  folding  sequence  in  order  to  form  a  stable  and  fully  packed  structure. 
Taking  a  USPS  mailbox  (Figure  2.2.1-7a)  as  an  example,  one  can  take  a  sheet  of  hard  board  paper  cut  into 
a  specific  shape.  One  folds  the  box  in  a  sequential  manner;  the  red  lines  represent  the  first  group  of  folds, 
where  some  lips  are  folded  to  form  supports;  the  blue  lines  represent  the  second  group  of  folds  to  partially 
close  the  box;  the  final  folds  (black  lines)  insert  the  two  small  lids  (on  the  left  top  and  bottom)  into  the 
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Figure  2.2. 1-5.  The  self  sequential  self-folding  (a)  without  collisions;  (b)  with  collisions.  When  there  is  no 
collision,  all  collision  indices  are  out  of  span  (0,l)(shown  in  (a)),  and  when  there  is  a  collision,  the  index 
falls  in  span  (0,  1)  (as  shown  in  (b),  two  points  (a  red  and  a  blue  point)  are  in  span  (0,1)).  The  red  dot 
between  (0,1)  corresponds  to  first  collision;  the  blue  dot  corresponds  to  the  second  collision. 


Figure  2.2. 1-6.  Schematic  of  the  interlocking  SMP  component  (a).  ROM  simulation  of  the  sequential 
folding  resulting  in  interlocking  (b),  the  collision  indices  of  the  chosen  design  (c),  and  a  comparison  of 
experiments  (plan  and  side  views)  and  the  ROM  simulations  (red  lines)  (d). 


holes  formed  by  the  long  blue  line  folds  to  increase  the  stability  of  the  box.  Inspired  by  this,  we  designed 
a  3D  folding  box  with  an  internal  locking  mechanism,  as  shown  in  Figure  2.2.1-7b.  Three  types  of  hinges 
were  used.  Hinge  H-3  is  used  on  the  smaller  lips  to  enable  fast  folding  of  these  lips  (red  lines).  Hinge  H-5 
is  used  in  the  location  indicated  by  the  green  line,  to  enable  intermediate  folding  speeds  such  that  some 
self-support  and  self-locking  mechanism  can  be  formed.  Hinge  H-6  (black  lines)  has  the  slowest  folding 
speed  and  is  used  to  form  the  final  structure.  The  original  shape  of  the  box  is  shown  in  Figure  2.2.1-7c. 
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Since  curved  portions  of  the  hinges  are  not  in  the  printing  plane,  these  hinges  are  suspect  to  damage 
during  unfolding.  To  improve  the  damage  tolerance,  hinges  are  specially  designed  (Figure  2.2.1-7d-f).  A 
periodic  step  profile  was  employed  along  the  surface  contour  of  the  hinge  (Figure  2.2.1-7d).  Thinner 
stepwise  sections  require  lower  bending  stress  and  have  higher  thermal  conductivity,  making  it  both 
easier  and  faster  to  unfold  complex  assemblies.  To  secure  the  interface  between  hinges  and  structural 
faces,  the  hinges  were  embedded  within  the  structural  faces  and  reinforced  with  a  triangular  interface  to 
maximize  axial  strength  (Figure  2.2.1-7e).  Longer  hinges  were  equipped  with  slots  cut  along  the  length  of 
the  jagged  contours  to  isolate  fracture  points  should  a  fracture  occur  (Figure  2.2.1-7f).  Figures  2.2.1-7g-j 
shows  the  folding  sequence.  The  whole  structure  folds  back  in  about  11  seconds. 

2.3.1  Multiple  shape  change  printed  active  composites 

In  this  work,  we  demonstrate  the  design  and  manufacture  of  active  composites  that  can  take  multiple 
shapes,  depending  on  the  environmental  temperature.  This  is  achieved  by  3D  printing  layered  composite 
structures  with  multiple  families  of  shape  memory  polymer  (SMP)  fibers  -  digital  SMPs  -  with  different 
glass  transition  temperatures  (Tg)  to  control  the  transformation  of  the  structure.  After  a  simple  single- 
step  thermomechanical  programming  process,  the  fiber  families  can  be  sequentially  activated  to  bend 
when  the  temperature  is  increased.  By  tuning  the  volume  fraction  of  the  fibers,  bending  deformation  can 


19 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


I  I  1  1  - 

44444444444444 


0 


j) 


Figure  2. 2. 1-7.  3D  folding  structures  mimicking  the  USPS  mailbox,  (a)  A  USPS  mailbox  is  folded  into  a  box 
by  following  a  sequence  of  folding;  b)  a  programmed  3D  printed  sheet  with  different  materials  assigned 
at  different  hinges;  c-f)  design  of  the  folding  box  with  some  details  at  the  hinges,  g-j)  upon  heating,  the 
sheet  folds  into  a  box  with  self-locking  mechanism. 


be  controlled.  We  also  develop  a  theoretical  model  to  predict  the  deformation  behavior  for  better 
understanding  the  phenomena  and  aiding  the  design.  This  work  is  in  the  manuscript  “Wu,  J .,  Yuan,  C., 
Ding,  Z.,  Isakov,  M.,  Mao,  Y.,  Wang,  T.,  Dunn,  M.L.,  Qi,  HJ.  2016.  Multi-shape  active  composites  by  3D 
printing  of  digital  shape  memory  polymers  Scientific  Reports  6,  Article  number:  24224" .  More  details  can 
also  be  found  in  this  manuscript. 


Design  concepts.  The  composites  in  this  study  consist  of  three  materials  with  different  glass  transition 
temperatures  (Tgs)  and  are  part  of  the  material  library  of  our  multi-material  3D  printer  (Objet260  Connex, 
Stratasys  Inc,  Edina,  MN,  USA).  Figure  2.3.1-la  shows  the  design  of  a  two-layer  composite.  The  matrix  is 
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Figure  2.3. 1-1.  Schematics  of  the  printed  SMP  composite  design,  (a)  The  design  of  the  two  layer  SMP 
composite  strips  and  the  characterization  of  the  design.  The  purple  color  represents  the  fiber  with  higher 
Tg  and  the  green  color  represents  the  fiber  with  lower  Tg.  (b)  The  typical  programming  steps  for  SMPs 
and  desired  response. 

TangoBlack+,  which  has  the  lowest  Tg  (~2°C)  of  the  three  materials.  Two  families  of  digital  SMP  fibers  with 
different  Tgs  are  embedded  in  the  two  layers,  respectively,  with  prescribed  volume  fractions.  The  fibers 
(fiber  1:  DM8530,  Tg~58°C;  fiber  2:  DM9895,  Tg~40°C,)  have  shape  memory  effects  in  the  temperature 
range  between  ~20°C  and  ~70°C.  The  thermomechanical  programming  steps  are  shown  in  Figure  2.3.1- 
lc.  We  first  stretch  the  composite  strip  at  an  elevated  programming  temperature  Th,  which  is  higher  than 
Tgs  of  both  fibers,  to  a  prescribed  strain  (e0),  and  then  cool  it  to  the  low  temperature  0°C  (TL,  which  is 
Lower  than  Tgs  of  both  fiber  materials)  while  maintaining  the  strain  e0.  The  applied  strain  is  then  released, 
which  fixes  the  composite  in  the  first  temporary  shape.  Upon  heating,  the  strip  gains  the  ability  to  deform 
due  to  the  temperature  dependent  viscoelastic  properties  of  the  matrix  and  fibers;  the  amount  of 
deformation  depends  on  the  thermomechanical  properties  of  the  fibers  and  the  matrix,  the  programming 
strain  and  the  ambient  temperature.  If  we  heat  the  sample  to  a  temperature  that  is  higher  than  the  Tg  of 
the  matrix  but  lower  than  that  of  the  fibers,  the  sample  is  able  to  change  to  a  second  temporary  shape. 
Due  to  the  long  stress  relaxation  time  of  the  fibers  at  low  temperatures,  we  are  able  to  obtain  a  series  of 
temporary  shapes  with  different  bending  curvatures  if  we  increase  the  temperature  in  a  staggered 
manner.  If  we  heat  the  sample  to  the  temperature  higher  than  the  Tgs  of  both  fibers,  it  will  recover  the 
flat  permanent  shape. 

The  shape  memory  behavior  of  the  SMP  composite  strip.  To  demonstrate  the  multi-shape  memory 
effects  of  the  3D  printed  SMP  composites,  we  print  the  active  composite  strip  by  following  the  design 
concept  illustrated  above.  Here,  the  strip  is  70mm  long,  6mm  wide,  and  2mm  thick,  with  the  fiber 
dimensions  of  a=b=0.38mm  (Figure  2.3.1-2a).  The  fibers  are  placed  at  the  center  of  each  layer,  i.e.,  the 
center  of  fiber  1  is  0.5mm  from  the  top  surface  and  the  center  of  fiber  2  is  0.5mm  from  the  bottom  surface. 
The  fibers  are  uniformly  spaced  at  0.5mm  (fiber  center  to  fiber  center)  to  give  volume  fractions  of  14% 
for  either  fiber  1  or  fiber  2.  In  order  to  facilitate  the  stretching  of  the  sample,  two  handles  are  directly 
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Figure  2. 3. 1-2.  Multi-shape  memory  effects  of  a  printed  active  composite  strip,  (a)  The  design  and 
dimensions  of  the  sample.  The  enlarged  drawing  is  the  cross  section  of  the  structure,  (b)  The  original 
printed  sample.  The  length  scale  in  the  bottom  is  in  mm.  (c)  -  (f)  Shape  change  of  the  sample  at  different 
temperature. 

printed  at  the  ends  of  the  sample,  and  one  handle  is  cut  off  before  reheating.  The  printed  sample  is  first 
stretched  at  high  temperature  TH=  70°C  with  strain  of  10%,  then  cooled  in  cold  water  of  ~2°C,  while  the 
loading  strain  is  kept.  Afterthe  sample  is  fully  cooled  down,  the  strain  loading  is  released.  Then  the  sample 

is  stepwise  heated  in  water  with  three  different  temperatures  (15°C,  30°C,  60°C)  to  demonstrate  the 
shape  recovery  process.  Figure  2.3.1-2b  shows  the  top  view  of  the  original  printed  sample.  Figures  2.3.1- 
2c-f  show  the  shape  memory  effects  of  the  sample.  Figure  2.3.1-2c  is  the  temporary  shape  of  the  sample 
at  2s  after  programming.  After  programming  the  sample  bends  slightly  towards  the  side  that  contains 
fibers  with  lower  Tg.  When  we  cool  for  one  minute  at  the  temperature  of  15°C,  which  is  lower  than  the  Tg 
of  fiber  2,  the  sample  bends  more  to  form  a  new  temporary  shape  (Figure  2.3.1-2d).  After  that,  we  heat 
the  sample  in  30°C  water  for  1  min. 

Figure  2.3.1-2e  shows  the  new  temporary  shape  of  the  sample.  Comparing  the  bending  angles  of 
the  sample  in  Figure  2.3.1-2d  and  Figure  2.3.1-2e,  we  can  see  the  shape  of  the  sample  changes  slightly 
although  the  temperature  is  higher  than  the  Tg  of  fiber  2.  Finally,  we  heat  the  strip  to  60°C  and  the  sample 
recovers  its  original  flat  shape  (Figure  2.3.1-2f).  The  sequence  of  transformations  of  the  3D  printed  SMP 
composites  indicates  that  we  can  achieve  multiple  shapes  using  convenient  3D  printing  technology  with 
a  programming  process  that  consists  of  a  single  simple  step.  In  addition.  Figure  2. 3. 1-2  shows  a  distinct 
feature  of  the  new  design:  unlike  conventional  SME  in  polymers  where  the  shape  change  is  simply  from 
one  to  another;  here,  a  simple  stretch  can  create  multiple  shape  changes,  i.e.  it  folds  (bends)  first,  then 
unfolds  (unbends)  to  become  straight. 

To  further  investigate  the  shape  recovery  process  and  speed  of  the  printed  active  composite  (PAC) 
strip,  we  directly  heat  the  sample  under  hot  water  to  prescribed  temperature  of  TH=70°C.  The  bending 
behaviorof  the  sample  can  be  described  by  the  variation  of  the  bending  curvature.  Figures  2.3. l-3a-f  show 
the  deformation  behavior  of  the  sample  at  different  times.  The  curvature  of  the  sample  as  a  function  of 
the  time  in  hot  water  is  plotted  in  Figure  2.3.1-3g.  It  is  clear  that  the  curvature  first  increases  when  the 
strip  begins  to  form  the  new  temporary  shape  and  then  decreases  when  the  sample  begins  to  recover  to 
its  permanent  shape  after  the  temperature  rises  higher  than  the  Tgs  of  both  fibers.  Note  that  the  whole 
recovery  time  is  only  12  seconds,  indicating  the  fast  response  speed  of  shape  change. 
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Figure  2.3. 1-3.  Bending  actuation  of  the  printed  active  composite  strip  during  the  recovery  process,  (a-f) 
Snapshots  of  the  sample  bending  at  different  times.  All  six  images  are  the  side  view  of  the  sample,  (g) 
Variation  of  the  bending  curvature  during  the  recovery  process. 


Theoretical  prediction  of  multi-shape  memory  behavior.  In  order  to  better  understand  the  deformation 
behavior  of  the  SMP  composites,  a  theoretical  model  is  developed  that  can  be  used  to  design  SMP 
composites.  In  this  model,  the  bending  behavior  of  the  composite  strip  is  described  using  classical 
laminate  theory,  and  a  multi-branch  model  is  used  to  describe  the  viscoelastic  mechanical  behavior  of  the 
matrix  and  fibers.  Since  the  recovery  speed  is  relatively  fast,  the  heat  transfer  problem  is  considered  to 
capture  the  non-uniform  temperature  distribution  in  the  sample  during  the  heating  process.  Using  the 
theoretical  model,  we  can  predict  the  deformation  behavior  of  the  SMP  composites.  The  curvature 
variation  of  the  SMP  composites  strip  is  calculated  as  a  function  of  the  heating  time.  Figure  2.3.1-4a  shows 
the  comparison  of  the  experimental  and  theoretical  results  of  the  curvature  variation  of  the  PAC  sample 
during  the  heating  process.  Overall,  the  model  captures  the  variation  of  the  sample's  curvature  during  the 
recovery  and  confirms  the  qualitative  explanation  of  the  operant  phenomena  of  the  preceding  section. 
The  model  predicts  a  recovery  that  is  about  1-2  seconds  faster;  this  might  be  because  the  strip  recovers 
in  water,  which  may  provide  some  resistance  to  the  recovery  motion. 

The  theoretical  model  provides  an  effective  tool  to  help  design  the  printed  active  composites.  For 
example,  the  programming  strain  influences  the  residual  strain  or  stress  stored  in  the  programmed  PACs. 
Using  the  same  dimensions  of  the  composite  shown  in  Figure  2.3.1-4a,  we  vary  the  prescribed  strain  and 
measure  the  curvatures  of  the  sample  in  the  recovery  procedure.  The  length  of  the  sample  L=70mm  and 
the  programming  temperature  TH=70°C,  TL=0°C.  Figure  2.3.1-4b  shows  the  initial  curvature  and  the 
maximum  curvature  of  the  PAC  as  a  function  of  the  pre-stretch  strain,  using  values  from  both  experiments 
and  theoretical  predictions.  It  is  clear  that  the  curvature  increases  with  the  pre-stretch  because  the 
residual  mismatch  strain,  which  is  the  driving  force  of  the  bending,  increases  with  the  programming  strain. 
The  effects  of  the  volume  fraction  of  the  fibers  are  investigated  by  experiments  and  simulations.  We  vary 
the  number  of  fibers  in  fiber  family  1  (which  has  the  higher  Tg  (in  purple  color))  to  vary  the  volume  fraction 
from  14.4%  to  7.22%,  3.61%  and  2.4%.  Comparison  between  the  experimental  results  and  theoretical 
predictions  of  the  initial  and  maximum  curvatures  are  plotted  as  a  function  of  the  volume  fraction  of  the 
fiber  with  higher  Tg  in  Figure  2.3.1-4c,  which  shows  that  when  the  volume  fraction  increases  the  maximum 
curvature  will  increase  monotonously;  however,  the  initial  curvature  is  nearly  unchanged. 


Printing  self-folding  and  self-opening  structures.  As  described  above,  by  tuning  the  volume  fraction  of 
fibers  and  the  prescribed  strain,  the  bending  of  the  SMP  composites  can  be  controlled.  Flere  we  show  an 
example  of  self-folding  and  opening  structures  that  exploit  this  phenomenon  and  our  understanding  of  it. 
More  examples  can  be  found  in  the  manuscript. 

We  first  design  a  self-assembling  and  disassembling  trestle.  Figure  2.3.1-5a  shows  the  design  of 
the  trestle,  which  has  4  identical  active  composite  strips  connected  at  the  center.  Figure  2.3.1-5b  shows 
the  cross-section  of  the  composite  strip.  Every  strip  has  one  fiber  with  higher  Tg  (fiber  1)  and  twelve  fibers 


23 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


5  10 

Volume  Fraction  (%) 


15 


Figure  2.3. 1-4.  Comparison  of  the  curvatures  between  experiments  and  theoretical  model,  (a)  Bending 
curvature  of  the  SMP  composite  vs  heating  time,  (b)  Effects  of  the  programming  strain  and  (c)  the  volume 
fraction  of  fibers  on  the  initial  and  maximum  curvature. 


with  lower  Tg  (fiber  2)  and  the  dimensions  of  the  strip  are  55mm(length)  by  6mm(width)  by 
2mm(thickness).  The  volume  fractions  of  fiber  1  and  fiber  2  are  1.2%  and  14.4%  respectively.  The 
structure  is  stretched  with  a  strain  of  8%  at  70°C  then  cooled  to  0°C.  After  releasing  at  the  low 
temperature,  the  sample  bends  slightly  with  a  small  curvature,  as  shown  in  Figure  2.3.1-5c.  When  heated, 
the  strips  fold  and  the  trestle  stands  up  from  the  flat  shape  as  shown  in  Figures  2.3.1-5d-f.  If  we  continue 
heating  the  sample,  the  trestle  goes  back  to  its  flat  shape  (Figure  2.3.1-5g).  This  design  can  be  used  as  an 
active  supporting  trestle  controlled  by  external  heating. 

Using  3D  printing  technology  to  print  shape  memory  polymers  provides  an  easy  and  efficient  way 
to  manufacture  active  structures,  and  of  course  it  can  be  used  to  create  complex  shapes  tailored  for 
specific  applications.  Flere  we  further  demonstrate  one  example:  a  smart  hook  (Figure  2.3. 1-6).  In  this 
design,  two  3D  printed  composites  strips  are  connected  at  the  ends.  After  being  stretched  with  10%  strain, 
cooled  to  0°C  and  relaxed  under  0°C  water  for  1  min,  one  end  handle  of  the  structure  is  cut  off.  Then  the 
nearly  flat  sample  can  be  used  as  a  hook  to  lift  a  small  box.  We  first  put  the  sample  into  hot  water  with 
temperature  of  30°C  and  the  straight  sample  bends  to  form  two  half  circular  shapes  (Figure  2.3.1-6b-c). 
Using  the  two  half  circular  shapes,  we  can  lift  up  a  small  basket  from  water  (Figure  2.3.1-6e).  To  release 
the  box  into  another  position,  we  increase  the  temperature  of  the  sample  to  be  higher  than  the  Tgs  of  the 
fibers  (Figure  2.3.1-6f-h).  From  these  simple  designs,  we  can  see  the  great  potential  of  using  the  3D  printed 
active  materials  in  creating  smart  structures  that  can  deform  on  demand  by  controlling  the  environment 
temperature. 
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Figure  2. 3. 1-5.  Self-assembling  and  disassembling  trestle,  (a)  The  design  of  the  trestle,  (b)  The  cross 
section  of  the  composites  strip,  (c)  The  shape  of  the  structure  after  programming,  (d)-(g)  The  deformation 
behavior  of  the  structure  in  the  recovery  process. 


Figure  2.3. 1-6.  Smart  hook.  (a)The  programmed  hook.  (b)-(c)The  bending  deformation  of  the  structure 
under  hot  water  with  temperature  of  30°C.  (d)-(e)  A  small  box  is  lifted  up  from  water.  (g)-(h)Releasing  the 
small  box  into  another  container.  The  water  in  the  container  is  in  70°C. 


2.4.1  Reversible  shape  changing  components  using  shape  memory  polymers  and 
hydrogels 

In  this  work,  we  demonstrate  a  new  reversible  shape-changing  component  design  concept  enabled  by  3D 
printing  two  stimuli-responsive  polymers  -  shape  memory  polymers  and  hydrogels  -  in  prescribed  3D 
architectures.  Our  approach  uses  the  swelling  of  a  hydrogel  as  the  driving  force  for  the  shape  change, 
and  the  temperature-dependent  modulus  of  a  shape  memory  polymer  to  regulate  the  time  dependence 
of  such  shape  change.  Controlling  the  temperature  and  aqueous  environment  allows  switching  between 
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Figure  2.4.1-1.  The  Schematic  graphs  of  the  design  concept  of  the  two-way  reversible  actuator,  a)  The 
reversible  actuation  component  by  3D  printing  where  the  hydrogel  is  confined  by  the  SMP  and  the 
elastomer  layers;  b)  The  key  dimensions  in  the  design,  c)  Schematic  plots  for  a  typical  two-way  actuation: 
In  SI,  the  sample  is  immersed  in  water  at  low  temperature;  in  S2,  it  is  brought  to  a  high  temperature 
environment  and  bends;  in  S3,  it  is  cooled  down  to  a  low  temperature;  in  S4,  it  is  allowed  to  dry  for  an 
elongated  time  period;  in  S5,  it  is  heated  to  recover  the  original  printed  shape.  This  finishes  one  actuation 
cycle,  which  can  be  repeated  many  times. 


two  stable  configurations  -  the  structures  are  relatively  stiff  and  can  carry  load  in  each  -  without  any 
mechanical  loading  and  unloading.  This  work  is  in  the  manuscript  "Mao,  Y.,  Ding ,  Z.,  Yuan,  C.,  Ai,  S.,  Isakov, 
M.,  Wu,  J.,  Wang,  T.,  Dunn,  Ml.,  and  Qi,  H.J.,  2016.  3D  Printed  Reversible  Shape  Changing  Component 
with  Stimuli  Responsive  Materials.  Scientific  Reports  6,  Article  number:  24761."  More  details  can  be  found 
in  this  manuscript. 

Design  concepts.  The  key  concepts  to  our  design  are  to  convert  the  hydrogel's  hydraulic  swelling  force 
from  equiaxial  to  a  linear  or  planar  force  that  can  drive  the  shape  change  in  one  particular  direction  or 
one  particular  plane,  and  to  use  the  temperature  sensitivity  of  the  SMP  properties  to  regulate  the  time 
for  such  shape  change.  Figure  2.4.1-la  shows  the  design,  where  hydrogel  and  elastomer  columns  are 
sandwiched  between  a  layer  of  the  SMP  (top)  and  a  layer  of  the  elastomer  (bottom).  Small  holes  are 
placed  in  the  elastomer  layer  to  allow  water  (or  other  solutions)  to  flow  in  and  out.  Figure  2.4.1-lb  shows 
the  notations  of  dimensions  in  the  design.  Figure  2.4.1-lc  shows  the  flow  chart  of  working  steps  for  a 
reversible  actuation  cycle.  The  printed  component  is  straight  after  printing.  It  is  then  immersed  into  water 
at  a  temperature  ~0°C  for  a  certain  amount  of  time  to  allow  the  hydrogel  to  absorb  the  water  (step  SI); 
in  addition,  due  to  the  low  temperature,  the  stiffness  of  the  SMP  is  high;  therefore,  the  volume  swelling 
of  the  hydrogel  is  highly  constrained,  and  the  strip  does  not  show  significant  shape  change.  Next  (step  S2), 
the  strip  is  brought  into  a  high  temperature  environment  (such  as  a  water  bath)  where  the  SMP  softens 
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significantly.  This  reduction  in  the  SMP  stiffness  allows  the  large  shape  change.  In  addition,  the  elastomer 
column  that  connects  top  and  bottom  layers  imposes  constraints  on  hydrogel  swelling  in  the  z-direction, 
thus  converting  the  swelling  force  into  the  x-y  plane.  Because  of  the  stiffness  difference  between  the 
elastomer  and  the  SMP,  the  strip  bends.  After  actuation  at  high  temperature,  we  cool  the  strip  to  a 
temperature  below  the  Tg  of  the  SMP  (step  S3);  due  to  the  increase  of  the  SMP  stiffness,  the  strip  is  stiff. 
In  the  ambient  environment,  the  hydrogel  loses  water  and  dries  in  step  S4.  After  the  hydrogel  is  fully  dry, 
immersing  the  strip  into  a  high  temperature  environment  recovers  the  straight  shape  of  the  strip  (step 
S5).  At  low  temperature,  the  strip  becomes  stiff  again.  This  finishes  one  cycle  of  reversible  actuation, 
which  can  be  repeated  multiple  times. 

We  used  a  multiple  material  3D  printer  (Objet260  Connex,  StrataSys,  Eden  Prairie,  MN,  USA)  to 
realize  the  above  concepts.  We  mainly  used  three  materials  in  the  printer  material  library:  Grey60, 
Tangoblack  (TB)  and  a  hydrogel.  As  shown  in  the  Method  Section,  Grey60  is  a  digital  material  with  the 
glass  transition  temperature  (Tg)  at  ~48°C,  which  can  be  used  as  the  SMP  when  the  temperature  is 
changed  between  0°C  and  60 °C.  TB  is  rubbery  at  the  RT  and  is  used  as  the  elastomer;  the  printed  hydrogel 
absorbs  water,  showing  a  RT  linear  swelling  ratio  of  1.18,  or  volumetric  swelling  ratio  of  1.64. 

Demonstration  of  the  two-way  actuation  of  a  simple  actuator.  The  concept  described  above  was 
implemented  into  a  simple  design,  which  was  then  printed  by  the  3D  printer.  Here,  the  top  SMP  layer 
was  0.5  mm  in  thickness,  and  the  bottom  TB  layer  was  1.5  mm,  and  hydrogel  layer  was  0.5mm.  The  radius 
of  the  hole  was  0.4mm.  The  pillar  had  a  square  cross-section  with  the  edge  length  of  0.5mm. 

The  printed  strip  is  shown  in  Figure  2.4.1-2a.  It  was  first  immersed  in  cold  water  (3°C)  for  12hrs 
and  showed  a  small  amount  of  bending  (Figure  2.4.1-2b).  The  strip  was  then  immersed  in  hot  water 
(temperature  of  75°C)  and  it  bent  into  the  shape  showing  in  Figure  2.4.1-2c  in  ~10  sec.  The  strip  was  then 
removed  from  the  water  and  was  cooled  down  to  the  RT  (Figure  2.4.1-2d).  The  strip  maintained  the 
bending  shape  (Figure  2.4.1-2e).  At  the  low  temperature,  it  was  stiff,  which  is  evident  in  Figure  2.4.1-2f 
where  the  bent  strip  can  support  a  deadweight  of  25g.  After  drying  at  low  temperature,  the  strip  was 
immersed  in  high  temperature  water  again  and  it  recovered  the  flat  shape. 


a)  The  printed  strip  is  straight;  b)  it  bends  slightly  after  immersing  in  cold  water  for  12hrs;  c)  it  bends 
quickly  when  immersing  in  hot  water;  d)  it  is  then  taken  out  of  water  to  the  room  temperature  air;  e)  it  is 
then  air  dried;  f)  the  strip  is  stiff  and  can  carry  a  load  of  25g.  If  it  is  heated,  it  returns  to  the  straight  shape. 


27 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Theoretical  studies  of  the  design.  Finite  element  method  (FEM)  simulations  were  conducted  to 
investigate  the  complicated  thermomechanical  coupling  of  the  SMP  and  the  swelling  property  of  the 
hydrogel  by  using  ABAQUS  (Dassault  Systems,  Johnston,  Rl,  USA).  The  thermomechanical  properties  of 
the  elastomer  and  the  SMP  were  considered  using  a  user  material  subroutine  based  on  the  finite  element 
constitutive  models  developed  by  the  authors.  The  non-equilibrium  swelling/shrinking  behaviors  of 
hydrogels  were  modeled  by  using  user  element  (UEL)  provided  by  Chester  et  al.  [11,  12].  To  reduce  the 
FEM  model  size,  only  a  quarter  of  a  unit  cell  was  considered. 

Figure  2.4.1-3b  shows  the  bending  angle  as  a  function  of  time  in  water.  It  is  clear  that  when  the  strip 
is  in  the  cold  water,  the  angle  evolves  slowly  as  water  is  absorbed  by  the  hydrogel.  However,  once  the 
strip  is  in  the  hot  water  (or  the  high  temperature  environment),  it  bends  quickly,  as  indicated  by  a  sudden 
jump  in  bending  angle  in  Figure  2.4.1-3b.  The  inset  in  Figure  2.4.1-3b  reveals  the  change  of  angle  occurs 
in  ~10s  in  hot  water. 
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Figure  2.4. 1-3.  Finite  element  modeling  of  water  in-take  and  actuation,  (a)  A  schematic  illustration  of  the 
finite  element  model,  (b)  Bending  angle  of  the  strip  as  a  function  of  time  of  the  strip  in  water.  The  inset 
shows  the  bending  angle  when  the  strip  is  immersed  in  hot  water.  Bending  angle  is  calculated  based  on  a 
strip  that  is  80mm  long,  (c)  FEA  simulation  results  of  the  weight  of  water  inside  the  hydrogel  as  a  function 
of  time.  Contour  plots  of  axial  stresses  in  the  matrix  and  hydrogel  and  the  swelling  ratio  in  hydrogels  at 
(d)  the  low  temperature  and  (e)  the  high  temperature. 
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Figure  2.4.1-3c  shows  the  water  absorbed  by  the  hydrogel.  It  also  shows  that  even  during  the  sudden 
change  of  bending  angle,  the  amount  of  water  in-take  only  changes  slightly,  indicating  that  the  modulus 
changes  in  the  SMP,  which  is  due  to  the  temperature  increase,  is  the  key  factor  for  the  increase  in  bending 
angles.  Figure  2.4.1-3d  shows  the  stress  contour  for  the  matrix  (left)  and  the  stress  and  swelling  ratio 
contours  for  the  matrix  (middle  and  right)  at  low  temperature  after  water  in-take.  It  can  be  seen  that  with 
the  water  in-take,  the  SMP  (top)  layer  develops  axial  stress  that  is  featured  in  bending  deformation 
(compression  on  the  top  and  bending  on  the  bottom).  In  addition,  the  swelling  of  hydrogel  is  non-uniform, 
with  most  swelling  occurring  in  the  area  where  hydrogel  is  in  direct  contact  with  water.  This  is 
understandable,  as  this  area  has  less  constraint.  Figure  2.4.1-3e  shows  the  same  contour  plots  at  the  high 
temperature.  In  comparison  with  the  contour  at  low  temperature,  there  are  two  differences.  First, 
although  the  trends  of  stress  distribution  in  the  SMP  are  the  same,  the  magnitude  decreases  significantly. 
This  is  because  of  the  significant  decrease  in  the  modulus  of  the  SMP.  Second,  the  volume  of  hydrogel 
increases.  This  is  interesting,  as  Figure  2.4.1-3b  shows  that  water  in-take  remains  almost  constant  as  we 
increase  the  temperature.  Therefore,  the  increase  of  hydrogel  volume  is  because  of  the  decrease  of 
hydraulic  pressure  and  the  compressible  nature  of  a  hydrogel.  It  should  be  noted  that  hydrogels  are 
usually  assumed  to  be  incompressible  because  their  shear  modulus  is  typically  much  smaller  than  their 
bulk  modulus.  This  assumption  is  true  when  the  hydraulic  pressure  is  not  high.  However,  under  the 
confined  environment  in  the  current  study,  the  significant  hydraulic  pressure  can  cause  observable 
volume  change  of  hydrogel.  To  further  confirm  that  the  sudden  bending  of  the  strip  is  caused  by  the 
change  of  the  modulus  of  the  SMP,  instead  of  water  in-take,  we  conducted  experiments  where  we  heated 
the  sample  by  using  an  oven.  The  results  show  that  although  the  bending  speed  is  slow,  the  strip  can  reach 
the  same  bending  angle  as  that  heated  in  hot  water. 

We  conducted  parametric  studies  on  the  effects  of  layer  thickness  to  the  bending  angle  and 
structure  stiffness.  The  details  about  these  studies  can  be  found  in  the  manuscript.  Here,  we  present  the 
case  of  parametric  study  on  actuation  speed.  Although  many  factors,  such  as  the  times  for  water 
absorption  and  depletion,  can  affect  the  actuation  speed,  we  demonstrated  the  cases  where  we  varied 
the  dimensions  of  the  component.  We  let  the  component  absorb  water  at  3°C  for  12  hours;  we  then  took 
it  to  the  high  temperature  of  75°C  to  allow  the  shape  change.  The  change  of  bending  curvature  as  a 
function  of  time  (measured  from  the  time  heating  starts)  during  heating  was  monitored  and  is  reported 
in  Figure  2.4. 1-4  for  the  components  with  different  layer  thicknesses,  represented  as  h  Elastomer*  h  hydrogel  •  hsMP 
in  the  figure.  Here,  we  varied  the  thickness  of  the  SMP  and  the  hydrogel  layers  but  maintained  the  total 
layer  thickness  to  be  constant  at  2.5mm  and  the  elastomer  layer  thickness  to  be  0.8mm.  It  can  be  seen 
that  a  thinner  SMP  layer  generally  promotes  a  larger  curvature  and  a  faster  rate  of  curvature  change.  The 


Figure  2.4.1-4.  Actuation  speeds  of  actuators  with  different  dimensions.  The  change  of  bending  curvature 
as  a  function  of  time  for  the  actuators  with  different  layer  thicknesses  (expressed  as  thicknesses  of 
Elastomer:Hydrogel:SMP)  during  the  Step  2  (heating). 
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difference  in  curvatures  and  curvature  change  rates  can  be  used  to  design  sequential  folding  patterns,  as 
will  be  shown  later  in  this  paper. 

2D  shape  change  objects.  To  demonstrate  the  potential  of  our  approach  on  reversible  actuation  of  the 
design,  we  used  the  bending  behaviors  to  create  different  shape  changing  structures.  In  the  first  example, 
we  applied  the  reversible  actuation  concept  to  a  hinge  design,  which  was  applied  to  a  ladder  where  some 
portions  of  its  structure  have  the  reversible  actuation  components  (Figure  2.4.1-5a).  The  hinge  of  the 
ladder  consisted  of  three  elements  discussed  above,  where  the  thicknesses  of  the  elastomer,  the 
hydrogel,  and  the  SMP  layers  were  0.5mm,  1.5mm  and  0.5mm  (shown  in  Figure  2.4.1-5b),  respectively. 
The  printed  ladder  is  shown  in  Figure  2.4.1-5c.  Similar  to  the  general  reversible  activation  cycle  discussed 
in  Figure  2.4.1-la,  after  allowing  the  ladder  in  the  low  temperature  water  (3  °C)  for  ~10  hrs,  the  hinges 
bent  slightly.  After  putting  it  in  high  temperature  water  (75  °C),  it  deformed  into  a  bench  as  shown  in 
Figure  2.4.1-5d.  The  bending  configuration  was  fixed  in  low  temperature  (3°C)  due  to  the  glass  transition 
property  of  the  SMP  layer.  After  drying  in  low  temperature  (3°C),  the  temporary  configuration  was 
sufficiently  stiff  and  could  sustain  a  50g  weight  as  shown  in  Figure  2.4.1-5e.  After  heating  the  bench  to 
75°C,  it  became  a  straight  ladder  again  (Figure  2.4.1-5f).  This  process  can  be  repeated  many  times. 

Reversible  folding/unfolding  origami.  Origami  folding  has  attracted  great  research  interest  in  recent 
years  due  to  its  potential  application  to  active  structures.  For  structural  engineers,  origami  has  proven  to 
be  a  rich  source  of  inspiration,  and  it  has  found  its  way  into  a  wide  range  of  structural  applications,  ranging 
from  wrapping  solar  cells  to  medical  stents  to  emergency  shelters.  Following  the  work  of  Fuchi  et  al.  [13], 
who  developed  a  design  methodology  for  origami  folding  actuators.,  an  origami  actuator  was  designed  by 
using  reversible  components.  Flere  the  hinge  element  had  the  thicknesses  of  0.25mm,  0.4mm  and 
0.15mm  for  the  elastomer,  the  hydrogel  and  the  SMP  layers,  respectively.  In  Figure  2.4.1-6a,  U  and  D 
indicate  the  bending  directions  of  the  hinge  with  U  meaning  upward  and  D  for  downward,  respectively. 
The  as-printed  sheet  is  shown  in  Figure  2.4.1-6c.  After  being  put  in  cold  water  for  lOhr  then  transferred 


Figure  2.4. 1-5.  A  reversible  ladder  and  bench  design,  a)  The  design  of  the  ladder  with  two  hinges;  b)  each 
hinge  contains  the  reversible  component,  c)  The  printed  ladder,  d)  It  bends  into  a  bench  shape  after 
swelling  in  low  temperature  water  for  lOhrs,  followed  by  heating  to  75°C  then  cooling  to  the  RT.  e)  The 
bench  is  stiff  and  can  carry  a  load  of  50g.  f)  The  bent  shape  is  recovered  when  the  bench  is  put  in  high 
temperature  water  (75°C).  The  scale  bars  represent  20mm. 
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Figure  2. 4.1-6.  Self-folding/unfolding  origami,  a)  A  schematic  graph  of  the  design;  b)  the  hinge  uses  the 
reversible  actuation  element;  c-e)  the  reversible  actuation  of  the  sheet  changing  from  flat  sheet  to  folded 
shape,  to  flat  sheet  again.  The  scale  bars  in  c-e)  are  45mm. 

to  hot  water,  the  plane  sheet  folded  into  a  3D  origami  structure.  After  being  sufficiently  dried  in  cold  air 
and  being  heated  by  hot  water,  the  sheet  recovered  its  initial  flat  configuration  (Figure  2.4.1-6e). 

As  a  final  example,  we  demonstrate  a  lotus  flower.  We  designed  and  printed  a  flower-shaped  3D 
structure  composed  of  three  types  of  petals,  and  experimentally  demonstrated  the  control  of  folding  of 
petals  as  a  function  of  controllable  geometrical  size.  According  to  the  discussion  in  Figure  2.4. 1-4,  for 
designs  with  the  same  total  and  elastomer  layer  thicknesses,  those  with  thinner  SMP  layers  promote 
larger  curvatures  and  faster  response  speeds.  By  elaborating  the  design  of  geometrical  size,  we  can  create 
the  sequence  folding  of  petals  like  a  flower.  The  three  layers  were  designed  with  different  thickness  ratios 
of  layers  of  the  elastomer,  the  hydrogel,  and  the  SMP  as  shown  in  Figure  2.4.1-7a.  The  inner  layer  was 
designed  with  the  thicknesses  of  the  elastomer,  the  hydrogel  and  the  SMP  layers  as  0.2mm,  0.4mm,  and 
0.3mm.  The  second  layer  was  designed  with  the  thicknesses  of  0.2mm,  0.3mm,  and  0.4mm  for  the 
elastomer,  the  hydrogel  and  the  SMP  layers;  and  the  outer  layer  was  designed  with  the  thicknesses  of  the 
elastomer,  hydrogel  and  SMP  layers  as  0.2mm,  0.2mm,  and  0.5mm.  The  corresponding  curvature  change 
profiles  were  characterized  by  three  different  speeds  for  three  layers  as  shown  in  Figure  2.4.1-7b.  After 
putting  this  structure  in  low  temperature  water  for  12hrs,  the  inner  layer  of  petals  bent  a  little  as  shown 
in  Figure  2.4.1-7c.  The  structure  was  then  immersed  in  high  temperature  water.  Immediately,  all  the  layers 
bent,  forming  a  flower-like  configuration  (Figure  2.4.1-7d  and  e).  Taking  the  structure  out  of  hot  water 
and  letting  it  dry,  the  structure  maintained  the  flower  shape  and  was  stiff.  As  shown  in  Figure  2.4.1-7g,  it 
could  carry  a  load  of  25g.  The  flower-like  structure  was  then  put  into  hot  water  and  the  structure  became 
flat  again  (Figure  2.4.1-7f).  This  process  can  be  repeated  many  times. 

2.5.1  Effects  of  oxygen  to  the  interface  during  the  layer-by-layer  photocuring  process 

Photopolymerization  is  one  of  the  most  widely  used  methods  for  additive  manufacturing  and 
microfabrication  of  polymer  structures.  However,  the  mechanical  properties  of  these  materials,  formed 
incrementally  or  layer-by-layer  by  photopolymerization,  remain  unclear.  One  critical  issue  is  the  strength 
of  the  interfaces  between  adjacent  layers.  During  free  radical  photopolymerization,  these  interfaces  are 
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Figure  2.4.1-7.  A  self-folding/unfolding  flower,  a)  Schematic  of  two  activated  shape  memory  petal-like 
structure;  b)  The  configuration  change  under  thermal  activation  as  functions  of  time;  c-f)  the  sequence  of 
reversible  actuation;  g)  the  dried  configuration  is  stiff  and  can  carry  a  load  of  25g.  The  scale  bar  in  g)  is 
12.5mm. 

exposed  to  atmospheric  oxygen,  which  is  detrimental  to  the  polymerization  reaction  due  to  radical 
inhibition.  The  influence  of  oxygen  on  the  interfacial  properties,  however,  is  still  not  well  understood.  In 
this  work,  the  role  of  oxygen  on  interfacial  strength  during  incremental  photopolymerization  is 
investigated.  We  carried  out  experiments  on  interfacial  bonding  under  different  photopolymerization 
processing  parameters  to  investigate  the  effect  of  oxygen  on  interfacial  strength.  Based  on  experimental 
results,  an  interfacial  model  was  developed  to  capture  the  evolution  of  interfacial  strength  under  oxygen 
inhibition.  This  work  was  published  in  "Zhao,  Z,  Mu,  X.,  Wu,  J.,  Qi,  H.J.,  Fang,  D.,  2016.  Effects  of  oxygen 
on  interfacial  strength  of  incremental  forming  of  materials  by  photopolymerization.  Extreme  Mechanics 
Letters  9(1):108-118." 

Methods:  The  resin  used  in  this  work  was  a  mixture  of  99.28  wt%  PEG-DA  (Polyethylene  glycol)  (700) 
diacrylate,  0.67  wt%  photoinitiator  Irgacure  819  (Phenylbis  (2,4,6-trimethylbenzoyl)phosphine  oxide  and 
0.05  wt%  photoabsorber  Sudan  I.  Realizing  the  challenge  in  conducting  mechanical  experiments  on  the 
single  interface  that  is  formed  in  the  layer-by-layer  photopolymerization  method,  we  used  an  alternative 
approach,  or  the  part-by-part  incremental  photopolymerization  method,  which  is  shown  schematically  in 
Figure  2.5. 1-1.  Samples  were  cured  between  two  glass  slides  with  1mm  thick  glass  supports  at  two  edges. 
At  first,  one  half  of  the  sample  (7.5mmxl0mmxlmm)  was  cured  for  different  time  periods  under  various 
light  intensities.  Before  photopolymerization  of  the  first  part,  an  oxygen  permeable 
Poly(dimethylsiloxane)  (PDMS)  strip  (cross  section  lmmx0.5mm)  (Figure  2.5.1-la)  or  oxygen  non- 
permeable  glass  (Figure  2.5.1-lb)  was  inserted  at  the  center  of  the  glass  mold.  These  two  cases  provide 
comparative  studies  on  the  effects  of  oxygen  to  the  interface  properties.  After  curing  the  first  part,  the 
strip  was  removed,  and  the  resin  was  injected  into  the  second  half  of  the  mold.  The  whole  structure  with 
a  total  length  of  15mm  was  then  cured  for  a  prescribed  sufficiently  long  time.  Tensile  tests  in  the  direction 
perpendicular  to  the  interface  were  carried  out  by  using  a  dynamic  mechanical  analysis  (DMA)  tester. 
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Figure  2.5. 1-1.  Part-by-part  photopolymerization  process,  (a)  Curing  of  the  first  part  by  inserting  a  PDMS 
lid.  (b)  Curing  of  the  first  part  by  inserting  a  glass  lid.  (c)  Curing  of  the  whole  structure  after  the  addition 
of  a  second  part. 

Fourier  transform  infrared  (FT-IR)  spectrum  of  the  sample  surface  was  taken  on  an  FT-IR  spectrometer 
with  an  attenuated  total  reflection  (ATR)  unit. 

SEM  and  microscopic  images  of  the  interfacial  structures  are  shown  in  Figure  2.5. 1-2.  In  these  figures, 
the  locations  of  the  interfaces  are  identified  by  the  ridges  and  are  marked  by  thin  dashed  lines.  It  is  clear 
there  are  no  apparent  discontinuities,  such  as  voids  or  cracks,  at  the  interfaces,  regardless  of  whether  the 
first  part  was  covered  by  PDMS  or  by  glass. 

Theoretical  model:  The  model  is  mainly  based  on  classical  methods  used  in  previous  work  on 
photochemistry  and  on  photo-mechanics,  but  several  additional  assumptions  are  introduced  to  further 
simplify  the  model.  The  modeling  domain  is  shown  in  Figure  2. 5. 1-3. 

Evolution  of  chemical  species  during  first  part  curing, 


Photoinitiator: 


photoinitiator  consumption  K - v - / 

photoinitiator  diffusion 


(2.5. 1-1) 


Free  Radical: 


dt 


(2.5. 1-2) 


radical  generation  self -termination  oxygen  inhibition 


Molecular  oxygen: 


o-PEGDA 


d2Co(*,0 

ax2 


(2. 5. 1-3) 


oxygen  consumption 


oxygen  diffusion 


Oxygen  diffusion  in  the  PDMS  strip. 
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Figure  2.5. 1-2.  SEM  and  optical  micrographs  of  the  interface,  (a)  SEM  image  of  sample  interface  prepared 
by  inserting  a  glass  strip  (b)  SEM  image  of  sample  prepared  by  inserting  a  PDMS  strip  (c)  Microscopic  image 
of  sample  interface  prepared  by  inserting  a  glass  strip  (d)  Microscopic  image  of  sample  interface  prepared 
by  inserting  a  PDMS  strip.  (Dashed  lines  were  added  to  indicate  the  position  of  interface). 


Figure  2.5. 1-3.  Schematic  graphs  of  the  modeling  domains  of  first  part  curing  (a)  and  second  part  curing 
(b). 
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Double  bonds  in  crosslinker: 


dCu(x,t) 

dt 


=  -kDCM(x,t)CR(x,t) 


pM 
double  bonds  consumption 

dco(x,t)=D  d2co(x,t) 


dt 


o-PDMS 


Evolution  of  chemical  species  during  second  part  curing, 


dr 

dCg(x,r) 

dr 


dxz 


=  2  J3C1  (x,  t)  -  ktC\  (x,  r)  -  k0CR  (x,  t)Cb  (x,  t)  , 


dC0(x,r)  =-k  C  (x  t)C  (x  t)  +  D  — 

~  L  )'^o\A/->  L  ^o-PEGDA  ~  2  ' 

or  ox 


dCM(x,r) 

dr 


=  ~kpCM(x,T)CR(x,T )  . 


(2.5. 1-4) 

(2.5. 1-5) 

(2.5.1-6a) 

(2.5.1-6b) 

(2.5.1-6c) 

(2.5.1-6d) 


Interfacial  bridging  reaction, 

=  kbCM  (0- ,  t)Cm  (0  +,t)  (2.5. 1-7) 

dr 

The  relation  between  interfacial  bonds  and  interfacial  strength  is  described  by  the  scaling  relation  crs  = 
aCb. 

Results:  Oxygen  effect  under  different  processing  conditions:  Interfacial  strength  as  a  function  of  first 
part  curing  time  is  shown  in  Figure  2.5.1-4a,  and  the  surface  FT-IR  spectrum  of  the  first  part  is  shown  in 
Figure  2.5.1-4b.  Two  observations  can  be  made  based  on  these  results.  First,  the  interfacial  strength 
decreases  with  the  curing  time  of  the  first  part,  because  there  are  fewer  unconverted  double  bonds  in  the 
first  part.  Second,  the  existence  of  oxygen  can  improve  the  interfacial  strength,  especially  for  longer  curing 
times.  This  comes  from  the  fact  that  oxygen  inhibition  can  decelerate  the  consumption  of  double  bonds 
during  first  part  curing.  Interfacial  strength  as  a  function  of  light  intensity  is  shown  in  Figure  2.5.1-4c. 
Although  the  presence  of  oxygen  improves  interfacial  behavior  even  if  light  intensity  reaches  15mWcm"2, 
the  interfacial  strength  shows  a  continuous  decrease  with  increases  in  light  intensity  as  a  result  of  higher 
reaction  rate.  Numerical  results  based  on  the  theoretical  model  are  shown  in  Figure  2.5.1-4a  and  Figure 
2.5.1-4c,  which  agree  well  with  experiments. 

Results:  Details  of  interfacial  reactions:  The  evolution  of  chemical  species  during  the  first  and  the  second 
curing  are  shown  in  Figure  2.5. 1-5.  Because  the  diffusivity  of  oxygen  in  PDMS  is  two  orders  of  magnitude 
higher  than  that  in  PEG-DA,  the  oxygen  concentration  at  the  resin  surface  (x  =  0  )  nearly  remains  the 
same  during  the  reaction  (Figure  2.5.1-5a  and  Figure  2.5.1-5b).  In  the  vicinity  of  this  surface  layer  where 
oxygen  inhibition  dominates  the  reaction  process,  the  number  of  reactive  radicals  is  significantly  reduced 
(Figure  2.5.1-5c).  As  a  result,  the  conversion  of  double  bonds  near  the  surface  is  reduced,  which  provides 
more  unreacted  sites  for  interfacial  bridging  upon  contact  with  the  second  part  (Figure  2.5.1-5d). 
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(a) 


(b)  (c) 

Figure  2.5. 1-4.  Interfacial  properties  of  the  two  parts  prepared  under  different  conditions,  (a)  Symbols: 
Tensile  strength  as  a  function  of  curing  times  for  the  first  part  (light  intensity  3.5mWcm"2);  Lines: 
concentration  of  interfacial  bridges  as  a  function  of  first  part  curing  time  (from  the  model  in  section  3). 
(b)  FT-IR  spectrum  of  liquid  PEG-DA  resin,  the  surface  of  the  first  part  covered  by  glass  or  PDMS  after  2min 
irradiation  and  after  lOmin  irradiation,  (c)  Symbols:  Tensile  strength  as  a  function  of  light  intensity  (curing 
time  of  first  part  2min);  Lines:  concentration  of  interfacial  bridges  as  a  function  of  light  intensity  (from  the 
model) 


The  interfacial  strength  shows  a  continuous  reduction  with  both  the  irradiation  time  and  the  light 
intensity,  as  a  result  of  a  decrease  in  the  number  of  double  bonds  in  the  first  part,  which  is  shown  in  Figure 
2.5.1-6a&b.  This  trend  is  decelerated  if  oxygen  is  present.  Because  radicals  are  significantly  reduced 
through  inhibition  reaction  (Fig.  2.5.1-5a),  many  more  double  bonds  are  left  unconverted,  which  is 
beneficial  to  interfacial  bridging  with  the  second  part.  In  Fig  2.5.1-6a,  for  the  case  without  oxygen, 
deviations  between  model  prediction  and  experiments  appear  at  longer  time  intervals.  This  may  be 
because  self-termination  no  longer  obeys  a  second  order  kinetic  rule  due  to  the  entanglement  effect  of 
the  network. 

In  order  to  better  illustrate  the  role  of  oxygen,  interfacial  bridging  is  calculated  by  varying  the  volume 
ratio  of  atmospheric  oxygen;  the  interfacial  strength  is  then  derived  by  using  the  scaling  relation 
<JS  —aCB.  As  can  be  seen  in  Figure  2.5.1-6c,  interfacial  strength  <j s  increases  by  almost  200%  as  the 

atmosphere  approaches  pure  oxygen.  The  strength  as  oxygen  volume  ratio  approaches  0  corresponds  to 
the  case  without  oxygen. 
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Figure  2. 5. 1-5.  Evolution  of  different  species  during  the  part-by-part  curing,  (a)  Oxygen  concentration  near 
the  interface  during  the  first  part  curing,  (b)  Oxygen  concentration  at  resin  surface  as  a  function  of  time 
and  PDMS  layer  thickness  (c)  Radical  concentration  near  the  free  surface  during  the  first  part  curing,  (d) 
Double  bond  concentration  in  two  sides  during  the  second  part  curing  (the  first  part  curing  time  is  5min). 
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Results:  Parametrical  studies:  We  then  conduct  parametrical  studies  by  varying  model  parameters.  Here, 
the  strengths  of  interfaces  with  ( <Jwslth )  and  without  oxygen  ( <j™lthouI )  are  calculated  respectively.  Figures 

2.5.1-7a  and  Figure  2.5.1-7b,  show  the  dependence  of  <Jwsuh  - cr‘v"/"""  on  resjn  composition,  specifically 
the  initial  photoinitiator  concentration  C/0  and  the  initial  double  bond  concentration  CM0.  The  red  dot 
in  Figure  2.5.1-7a  indicates  the  experimental  condition.  As  shown  in  both  Figures  2.5.1-7a  &b,  the  effect 
of  oxygen  is  more  apparent  with  the  increase  of  CM0  at  a  fixed  C/0 .  This  trend  can  also  be  found  in 
Figure  2.5.1-7c.  According  to  our  model,  the  conversion  rate  of  double  bonds  is  proportional  to  a  product 
of  CM  and  CR,  which  relate  to  CM0  and  C/0,  respectively.  For  a  constant  initial  photoinitiator 

concentration,  increasing  the  number  of  double  bonds  implies  that  more  double  bonds  will  be  left 
unreacted.  This  is  more  pronounced  in  the  presence  of  oxygen  at  the  interface,  as  oxygen  consumes  a 
significant  number  of  radicals.  This  analysis  can  also  be  used  to  explain  Figure  2.5.1-7d  where  the  initial 
photoinitiator  concentration  increases  while  maintaining  the  double  bond  concentration  constant.  In  this 
case,  the  increase  of  photoinitiators  generates  more  radicals,  which  consume  more  double  bonds  at  the 
interface  and  thus  reduce  the  interfacial  strength. 
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(a)  (b) 


(c) 

Figure  2. 5. 1-6.  (a)  Double  bond  concentration  near  the  surface  as  a  function  of  first  part  curing  time 
(incident  light  intensity  I0  =  3.5 mWcm  2 ).  (b)  Double  bond  concentration  near  the  surface  as  a  function 
of  incident  light  intensity  (first  part  curing  time  2min).  (c)  Interfacial  strength  as  a  function  of  atmospheric 
oxygen  volume  ratio  (first  part  cured  for  2min,  I0  =  3.5 mWcm~2 ). 


Figure  2.5. 1-7  suggests  that  a  resin  composed  of  a  relatively  small  number  of  photoinitiators  and  a  large 
number  of  double  bonds  (which  can  be  realized  by  using  short  chain  crosslinkers)  seems  desirable  for  the 
interfacial  strength  of  incrementally  photopolymerized  structures.  The  discussion  of  Figure  2. 5. 1-7  shows 
that  the  conversion  rate  of  double  bonds  plays  an  important  role  in  determining  the  interfacial  strength. 
In  general,  a  slow  conversion  rate  of  double  bonds  at  the  interface  in  the  first  part  curing  will  lead  to  a 
stronger  interface.  However,  a  slow  double  bond  conversion  rate  in  the  bulk,  which  can  be  controlled  by 
the  photoinitiator  concentration  and  the  initial  double  concentration,  is  not  desirable  in  layer-by-layer 
applications,  as  this  leads  a  condition  where  the  whole  material/structure  fails  to  fully  cure.  Oxygen,  which 
only  dominates  at  the  interface,  provides  an  efficient  method  of  reducing  the  double  bond  conversion 
rate  at  the  interface,  thus  it  is  beneficial  to  the  interfacial  strength  and  to  the  overall  properties  of  the 
part  or  structure. 

2.5.2  Origami  Folding  by  Frontal  Photopolymerization 

Origami  structures  are  of  great  interest  in  microelectronics,  soft  actuators,  mechanical  metamaterials  and 
biomedical  devices.  Current  methods  of  fabricating  origami  structures  still  have  several  limitations,  such 
as  complex  material  systems  or  tedious  processing  steps.  In  this  work,  we  present  a  simple  approach  for 
creating  3D  folding  origami  structures  by  the  frontal  photopolymerization  method,  which  can  be  easily 
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Figure  2.5. 1-7.  (a)  Strength  increase  due  to  oxygen  as  a  function  of  initiator  concentration  C/0  and  double 
bond  concentration  of  CM0  (first  part  cured  for  2min,  incident  light  intensity  3.5  mW  cm  2),  (b)  Strength 
increase  due  to  oxygen  as  a  function  of  CM0  /  C/0  (extracted  from  (a)),  (c)  Interfacial  strength  as  a  function 
of  double  bond  concentration  CM0  ( C10  =  0.02  mol  L~l ,  extracted  from  (a)),  (d)  Interfacial  strength  as  a 
function  of  initiator  concentration  C/0  ( CM0  =  6  mol  L~l ,  extracted  from  (a)). 


implemented  by  using  a  commercial  LED  projector.  This  work  is  published  in  "Zhao,  Z,  Wu,  J.,  Mu,  X.,  Qi, 
H.J.,  Fang,  D.,  2016.  Origami  Folding  by  Frontal  Photopolymerization,  submitted/' 

The  concept  of  our  method  is  based  on  the  volume  shrinkage  during  photopolymerization.  When 
a  polymer  film  is  cured  from  one  side  in  the  frontal  photopolymerization  process,  the  newly  cured  layers 
shrink  under  the  confinement  of  older  layers.  By  adding  photoabsorber  into  the  polymer  resin,  an 
attenuated  light  field  is  created  and  leads  to  a  nonuniform  stress  field  in  the  thickness  direction,  which 
drives  the  film  to  bend  towards  the  newly  formed  side.  The  degree  of  nonuniformity  can  be  controlled  by 
adjusting  the  LED  light  gray-scale  and  the  irradiation  time,  an  easy  approach  for  creating  origami 
structures.  The  behavior  is  examined  both  experimentally  and  theoretically.  Two  methods  are  also 
proposed  to  create  different  types  of  3D  origami  structures. 

The  concept  of  our  method  is  based  on  photopolymerization  induced  volume  shrinkage.  During 
photopolymerization  of  the  resin,  the  material  volume  decreases  as  a  result  of  covalent  bond  formation 
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Figure  2. 5. 2-1.  Shrinkage-induced  bending  during  frontal  photopolymerization.  (a)  Schematic  process  of 
volume-shrinkage-induced  bending:  Sequential  shrinkage  occurs  during  the  frontal  photopolymerization 
of  a  polymer  sheet,  and  internal  stress  developed  during  the  process  drives  the  sheet  to  bend,  (b) 
Schematic  process  of  shrinkage-induced  bending  test  (with  PDMS  substrate),  (c)  Schematic  process  of 
shrinkage-induced  bending  test  (without  PDMS  substrate). 


between  monomers  and  crosslinkers.  This  behavior  is  commonly  assumed  to  be  undesirable  in  additive 
manufacturing  (or  3D  printing),  microfabrication  and  manufacturing  of  polymer  composites,  because  it 
causes  shape  distortion  and  internal  stress.  In  frontal  photopolymerization,  the  moving  front  of  cured 
region  leads  to  a  nonuniform  volume  shrinkage,  which  can  be  controlled  and  utilized  to  create  folding 
structures  and  3D  origami  shapes 

The  process  is  shown  schematically  in  Figure  2.5.2-la.  In  order  to  control  light  penetration, 
photoabsorber  is  added  to  induce  an  intensity  gradient  in  the  resin.  Upon  irradiation,  the  surface  layer 
directly  exposed  to  light  will  be  cured  immediately,  while  the  resin  below  this  layer,  without  sufficient 
illumination,  still  remains  in  the  liquid  state.  As  the  front  of  liquid-solid  transition  propagates  toward  the 
liquid  resin,  the  thickness  of  the  cured  polymer  film  increases.  Due  to  the  presence  of  a  light  intensity 
gradient  caused  by  the  photoabsorber,  volume  shrinkage  happens  in  a  sequential  manner.  That  is,  early 
cured  materials  shrink  first,  while  subsequently  cured  parts  shrink  under  the  confinement  of  the  older 
parts.  As  a  result,  a  nonuniform  stress  field  develops  across  the  sample's  thickness,  with  a  tendency  to 
bend  towards  the  newly  formed  parts.  Because  of  the  restriction  imposed  by  the  uncured  liquid  resin,  film 
that  has  an  internal  stress  gradient  remains  flat  in  shape.  Bending  occurs  once  the  film  is  moved  out  of 
the  liquid,  and  is  related  to  several  processing  conditions,  such  as  light  intensity  and  illumination  time.  By 
imposing  spatial  control  over  the  curing  condition,  complex  3D  origami  structures  can  be  created. 

Results  and  discussions  -  Experimental  investigation  of  shrinkage  induced  bending:_The  property  of 
acrylate  resin  was  characterized  by  checking  photopolymerized  thin  films  with  different  illumination  times 

under  a  constant  intensity  of  SmWcnf2 ,  and  the  results  showed  the  Young's  modulus  of  the  sample 
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increases  linearly  with  the  normalized  conversion  degree  (p  of  C=C  double  bonds.  Significant  shrinkage 
strain  around  -5%  occurs  after  solidification,  and  its  variation  becomes  much  more  mild  as  the  reaction 
proceeds.  During  frontal  photopolymerization,  the  evolution  of  these  two  properties  can  be  combined  to 
control  the  development  of  the  internal  stress  field. 

In  order  to  demonstrate  the  concept  of  shrinkage  induced  bending,  polymer  strips  were  cured  in 
a  glass  petri  dish  filled  with  liquid  resin,  with  LED  light  projected  upwards  from  the  bottom.  Inspired  by 
the  lubrication  effect  from  both  residual  oxygen  inhibition  and  pendent  oligomers  at  the  surface  of 
poly(dimethylsiloxane)  (PDMS),  a  PDMS  substrate  was  inserted  to  make  sure  that  cured  parts  would  not 
stick  to  the  bottom  so  that  the  sample  could  be  removed  easily.  The  process  is  shown  schematically  in 
Figure  2.5.2-lb.  After  different  periods  of  irradiation  under  various  light  intensities,  the  solidified  strip  was 
moved  out  of  the  liquid  to  measure  its  bending  curvature.  Figure  2.5.2-2a  shows  some  selected  samples 

that  were  cured  under  light  intensity  of  5mWcm~2 .  The  rectangular  strips  were  cut  short  to  facilitate 
measurement.  The  bending  curvatures  were  measured  from  several  figures  similar  to  Figure  2.5.2-2a  and 
are  summarized  in  Figure  2.5.2-2b  in  symbols.  It  is  interesting  to  note  that  curvature  decreases  with  both 
illumination  time  and  incident  light  intensity.  Although  the  sequential  shrinkage  behavior  remains  the 
same,  sample  stiffness  increases  due  to  the  increased  thickness  and  modulus,  making  it  more  difficult  to 
bend.  Bending  almost  disappears  when  the  curing  time  is  longer  than  10s  under  the  intensity  of 
9.5 mWcnf2  . 

For  a  comparison  group,  the  PDMS  substrate  in  the  petri  dish  was  removed,  and  then  rectangular 

strips  were  cured  in  the  same  manner  with  incident  intensity  of  SmWcnf 2  (  Figure  2. 5. 2-1).  The  results 
are  shown  by  red  open  circles  in  Figure  2.5.2-2b.  In  this  case,  the  bottom  surface,  which  cured  earliest, 
stuck  to  the  substrate  and  could  not  shrink  freely,  as  is  the  case  with  a  PDMS  layer.  After  the  whole  part 
was  cured  and  peeled  off  from  the  substrate,  the  bottom  surface  shrank  to  release  internal  stress.  This 
tendency  is  in  contrast  with  the  fact  that  the  bottom  part  was  compressed  during  the  sequential  shrinkage 
process,  and  the  overall  bending  curvature  was  reduced.  Therefore,  the  boundary  condition  at  the  bottom 
of  the  substrate  plays  a  key  role  in  determining  shape  changing  behaviors.  The  effect  of  the  surface  layer's 


Figure  2. 5. 2-2.  Experimental  results  of  shrinkage-induced  bending,  (a)  Rectangular  samples  cured  in  a  petri 

dish  with  different  irradiation  time  (light  intensity  5 mWcnf2).  (b)  Bending  curvature  of 
photopolymerized  strips  as  a  function  of  irradiation  time  and  incident  light  intensity,  (c)  Stress  distribution 

across  sample  thickness  (light  intensity  5 mWcnf2 ). 
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confinement  becomes  negligible  when  the  total  thickness  of  cured  sample  increases,  because  the  overall 
stress  gradient  along  thickness  governs  the  bending  behavior. 

Theoretical  model  of  shrinkage  induced  bending:_The  shrinkage-induced  bending  is  described  by  a 
theoretical  model  based  on  beam  theory.  During  photopolymerization,  samples  are  constrained  by  liquid 
resin.  The  out-plane  rotation  is  restricted,  and  only  in-plane  stretching  occurs.  The  mechanical  behavior 
of  films  under  illumination  should  be  governed  by  force  equilibrium  along  the  thickness  direction, 

rz(t)  pz{t)  c 

crdz  =  E(s-ss)dz  =  0.  (2. 5. 2-1) 

Here  z  is  the  solidified  thickness.  The  stress  a  is  decided  by  <J  =  Ese  =  E(8  —  8S),  where  E  is  the 
Young's  modulus,  se ,  8  and  8s  are  the  elastic  strain,  the  total  strain  and  the  shrinkage  strain, 
respectively.  E  and  8s  are  dependent  on  the  conversion  degree  <p  of  double  bonds  in  the  acrylate  resin 
(details  of  the  photopolymerization  model  can  be  found  in  the  Supporting  information).  The  total  strain 
8  and  shrinkage  strain  8S  are  set  to  zero  if  the  conversion  degree  (p  is  below  the  gel  point  cpc .  It  should 

be  noted  that  the  integral  domain  of  Equation  (2. 5. 2-1)  varies  with  time.  Thus  it  is  solved  incrementally 
by  time  intervals  to  get  the  value  of  8  at  different  z  and  t. 

When  the  solidified  film  is  moved  out  of  the  liquid  resin,  the  rotation  constraint  disappears,  and 
the  part  will  bend  under  its  internal  stress  field.  Both  force  equilibrium  and  moment  equilibrium  should 
be  satisfied, 


pz(  t)  pz(t) 

|  adz=\  Eeedz  =  0 

Jo  Jo 

(2.5.2-2a) 

rz(t )  fz(t) 

f  zodz  =  \  zEs  dz  =  0 

Jo  Jo 

(2.5.2-2b) 

Here  the  elastic  strain  is  obtained  by  a  superposition  of  bending  and  stretching, 

Se  =£-£S  +£b+K(z-^) 

(2. 5. 2-3) 

where  8b  is  the  bending  strain  of  the  midplane  at  z  —  zl  2 ,  and  k  is  the  bending  curvature.  The  value 
of  8  for  in-plane  stretching  is  passed  in  from  the  calculation  of  Equation  (2. 5. 2-1).  Equation  (2.5.2-2a)  and 
(2.5.2-2b)  can  be  solved  as  a  whole  to  get  the  value  of  8b  and  k  . 

Theoretical  calculation  results  based  on  the  model  are  shown  in  Figure  2.5.2-2b  in  lines.  The 
dependence  of  the  curvature  k  on  the  irradiation  time  and  intensity  is  captured  well.  An  interesting 
finding  is  that,  as  a  result  of  more  uniform  intensity  distribution  and  higher  rate  of  thickness  increase,  the 
decrease  of  curvature  with  irradiation  time  is  faster  under  a  stronger  incident  light.  This  observation  leads 
to  a  design  concept  wherein,  by  varying  light  intensity  (or  gray  scale  pattern)  in  a  2D  illumination  pattern, 
we  are  able  to  create  thick  stiff  panels  or  soft  thin  bending  parts:  for  the  same  light  irradiation  time,  parts 
under  high  light  intensity  have  high  stiffness  but  no  bending,  and  parts  under  low  light  intensity  are 
relatively  soft  but  with  a  sufficient  bending  curvature. 

Based  on  Equation  (2. 5. 2-1)  in  the  model,  the  evolution  of  the  stress  field  for  a  sample  cured 

under  irradiation  of  5 mWcnf2  is  shown  in  2.5.2-2c.  The  bottom  layer  that  was  cured  at  first  stayed  in  a 
compressive  state,  while  the  top  part  was  under  tension.  The  absolute  value  of  stress  was  higher  in  the 
bottom  part  because  of  much  stiffer  material.  As  the  irradiation  time  increases,  although  a  stress  gradient 
still  exists,  the  thickness  and  modulus  of  the  cured  sample  becomes  sufficiently  large  to  resist  bending 
deformation,  and  the  overall  curvature  is  negligible. 
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Origami  structures  created  by  one-side  illumination^  Utilizing  shrinkage-induced  bending  in 
photopolymerization,  we  can  create  complex  3D  structures  from  flat  polymer  sheets  with  programmed 
2D  light  patterns.  It  should  be  pointed  out,  if  the  polymer  sheet  is  cured  by  illumination  from  the  bottom, 
as  we  explained  before,  it  can  only  bend  towards  the  top  side.  We  will  show  first  how  to  produce  origami 
structures  from  this  type  of  one-side  illumination.  Using  an  LED  projector,  2D  light  patterns  with  spatial 
control  over  intensity  are  created  through  variation  of  the  gray  scale.  Low  intensity  is  used  at  those  points 
where  folding  (or  large  bending)  is  expected.  For  the  table  structure  in  Figure  2.5.2-3a-d,  intensities  were 

set  to  15 mWcnf2  and  3 mWcnf2  for  the  panel  and  legs,  respectively.  This  was  realized  by  creating 
figures  with  a  spatial  variation  of  grayscale  (inset  in  Figure  2.5.2-3a).  After  a  specific  irradiation  time  (10s 
for  Figure  2.5.2-3a),  as  discussed  before,  the  cured  sheet  remained  flat  while  immersed  in  the  liquid  resin. 
Parts  that  were  exposed  to  high  intensity  were  thicker,  and  parts  exposed  to  low  intensity  (for  example 
those  legs  in  Figure  2.5.2-3a)  were  relatively  thinner.  When  the  sample  was  removed  from  the  liquid, 
those  parts  exposed  to  low  intensity  (legs  in  Figure  2.5.2-3b)  bent  towards  the  newly  formed  side  as  a 
result  of  nonuniform  stress.  Because  the  bending  parts  had  lower  stiffness,  this  configuration  was  not 
stable.  To  fix  the  shape,  samples  were  held  by  a  tweezer  and  further  cured  under  a  uniform  LED  light  of 

10  mWcnf2  for  20-30s  (Figure  2.5.2-3c),  which  stiffened  the  structure  due  to  the  dramatic  increase  of  the 
modulus  in  the  previously  soft  parts.  After  this  post  curing,  the  sample  became  strong  enough  to  support 
some  weight.  A  table  structure  holding  several  glass  cover  slips  (with  a  total  weight  of  l.lg)  is  shown  in 
Figure  2.5.2-3d. 


Figure  2. 5. 2-3.  Fabrication  of  origami  structures  by  one-side  illumination,  (a)  Polymer  sheet  right  after 
photopolymerization.  (b)  Free  bending  of  spatial  differently  cured  sheet,  (c)  Shape  fixing  of  bending 
structures  by  post  curing,  (d)  Sample  after  post  curing  is  stiff  and  able  to  hold  several  glass  cover  slips,  (e) 
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Flower  structures  with  different  opening  degree,  (f)  Polymer  sheet  with  a  continuous  variation  of 
curvature. 


a) 


c) 


ATA 


Figure  2. 5. 2-4.  Polyhedrons  created  by  one-side  illumination  (a)  Octahedron  (panel  intensity  15  mWcnf2 
,  hinge  intensity  3 mWcnf2 ,  irradiation  time  6s).  (b)  The  grayscale  pattern  of  octahedron  structure,  (c) 
Truncated  cube  (panel  intensity  15 mWcnf2 ,  hinge  intensity  3.3 mWcm  2,  irradiation  time  6s).  (d) 
Regular  icosahedron  (panel  intensity  15  mWcm~2 ,  hinge  intensity  4 mWcm~2 ,  irradiation  time  6s). 


Much  more  delicate  control  over  the  light  pattern  produces  structures  with  higher  complexity. 
Flowers  with  different  opening  degrees  are  shown  in  Figure  2.5.2-3e.  These  four  flowers  were  all  created 

using  an  irradiation  period  of  6s.  The  illumination  intensity  of  the  central  round  panel  was  10 mWcm  2, 

and  petal  intensities  were  set  to  2.5,  4,  6  and  7  mWcm  2 ,  respectively.  The  sample  that  was  illuminated 
with  the  lowest  intensity  shows  the  smallest  opening  angle.  The  round  panel  at  center  is  illuminated  with 
a  high  intensity,  and  it  remains  flat  in  all  4  cases.  A  polymer  sheet  with  continuous  variation  of  bending 
curvature  is  shown  in  Figure  2.5.2-3f.  This  structure  was  realized  by  illuminating  a  trapezoidal  pattern  with 
a  gradient  of  grayscale  for  4s.  The  edge  with  the  highest  curvature  was  illuminated  with  the  lowest 

intensity  (2.5  mWcm~2 ),  and  intensity  gradually  increased  till  the  other  edge  (to  10  mWcm~2). 

3D  folding  structures  were  created  by  using  a  combination  of  flat  panels  and  hinges  bent  towards 
one  side.  Some  different  types  of  polyhedron  capsules  are  shown  in  Figure  2. 5. 2-4.  In  order  to  stabilize 
the  folding  process,  extra  panels  are  introduced  in  the  light  pattern.  For  example,  to  form  the  shape  in 
Figure  2.5.2-4a,  four  more  panels  marked  by  red  crosses  are  attached  to  the  net  pattern  of  octahedron  in 
Figure  2.5.2-4b.  During  the  folding  process,  an  extra  panel  overlaps  on  the  top  of  a  neighbor  panel  and 
thus  stabilizes  the  structure. 


Origami  structures  created  by  two-sides  illumination:  Structures  presented  in  the  previous  section  are 
created  by  polymer  sheets  bending  towards  one  single  side,  thus  the  freedom  of  design  is  limited.  In  this 
section,  we  will  show  a  modified  method  that  is  able  to  create  more  complex  origami  structures.  The 
process  is  shown  schematically  in  Figure  2.5.2-5a.  Liquid  resin  is  injected  into  a  flat  mold  confined  by  two 
glass  slides,  two  PDMS  lubrication  layers  are  pasted  on  the  inner  sides  of  the  glass  slides  respectively  to 
avoid  adhesion  of  cured  parts.  The  distance  between  two  PDMS  slides  is  confined  to  0.5mm.  For  the  first 
step,  parts  designed  to  bend  towards  direction  1  (indicated  by  number  1  in  Figure  2.5.2-5a)  are  illuminated 
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under  low  intensity  from  the  bottom.  After  that,  the  whole  structure  is  flipped  over  without  in-plane 
movement  of  position.  Then  parts  designed  to  bend  towards  direction  2  (indicated  by  number  2  in  Figure 
2.5.2-5a)  are  illuminated  in  the  same  manner.  By  releasing  the  mold,  the  cured  sheet  bends  towards  the 
two  different  sides.  In  a  real  application,  one  can  also  apply  two  light  sources,  or  use  a  light  splitter  to 
expedite  the  processing  time.  Two  light  patterns  are  projected  in  sequence  to  realize  this  two-step 


a) 


ft 


Glass 
PDMS 
Liquid  resin 


ttttt  ft 

1st  light  pattern 


tftft  fttft 


2nd  light  pattern 


Figure  2. 5. 2-5.  Fabrication  of  origami  structures  by  two-sides  illuminations,  (a)  Schematic  process  of  two- 
sides  illumination  method,  (b)  Table  (panel  intensity  18  mWcnf2 ,  petal  intensity  2.5  mWcnf2 ,  irradiation 
time  of  each  side  7s),  (c)  Miura  structure  (panel  intensity  20 mWcnf2 ,  hinge  intensity  2 mWcnf2 , 
irradiation  time  of  each  side  7s),  (d)  tent  (panel  intensity  18 mWcm  2 ,  hinge  intensity  2  mWcnf2 , 
irradiation  time  of  each  side  7s)  and  (e)  crane  (panel  intensity  20 mWcnf2 ,  body  hinge  intensity  1 

mWcnf2 ,  wing  hinge  intensity  ImWcmf2 ,  irradiation  time  of  each  side  15s).  created  by  two-sides 
illumination  method. 
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illumination.  The  positions  marked  by  2  are  made  black  during  the  first  step  illumination,  and  the  positions 
marked  by  1  are  made  black  during  the  second  illumination,  while  plane  panels  are  irradiated  in  both 
steps. 

3D  origami  structures  produced  by  two-side  illumination  are  presented  in  Figure  2.5.2-5b-e. 
Figure  2.5.2-5b  shows  a  table  with  legs  in  two  directions,  created  by  illuminating  the  leg  parts  from  two 
sides.  The  intersection  points  of  different  hinges  are  left  as  voids  to  facilitate  the  folding  process.  This  is 
shown,  for  example,  by  the  vertices  in  the  Miura  structure  (Figure  2.5.2-5c)  and  the  top  of  a  pyramid  tent 
(Figure  2.5.2-5d).  Different  intensities  can  be  used  for  separate  hinges  to  enrich  the  freedom  of  the 

origami.  For  the  origami  crane  in  Figure  2.5.2-5e,  an  intensity  of  ImWcnf2  was  used  for  the  wing  hinges 

and  a  very  low  intensity  of  1  mWcm  2  was  used  for  the  hinges  of  legs  and  tails. 

A  comparison  between  these  two  methods  shows  that  they  are  suitable  for  different  conditions. 
While  the  second  method  of  two-side  illumination  has  the  ability  to  create  very  complex  structures,  extra 
processing  steps  become  inevitable,  such  as,  for  example,  the  preparation  of  a  mold  and  the  insertion  of 
liquid  resin.  As  for  the  one-side  method,  structures  can  be  created  by  just  shining  a  single  pattern  onto  a 
liquid  resin  vat,  but  the  degree  of  freedom  of  the  origami  is  confined.  Thus,  if  a  structure  can  be  designed 
simply  by  a  combination  of  hinges  turned  towards  one  direction,  the  one-side  method  seems  sufficient. 
The  two-side  illumination  method  becomes  necessary  when  more  delicate  control  over  the  origami  is 
desired. 

2.6.1  Thermal  cure  effects  on  electromechanical  properties  of  conductive  wires  by  direct 
ink  write  for  4D  printing  and  soft  machines 

Recent  developments  in  soft  materials  and  3D  printing  are  promoting  the  rapid  development  of  novel 
technologies  and  concepts,  such  as  4D  printing  and  soft  machines,  that  in  turn  require  new  methods  for 
fabricating  conductive  materials.  Despite  the  ubiquity  of  silver  nanoparticles  (NPs)  in  the  conducting 
electrodes  of  printed  electronic  devices,  their  potential  use  in  stretchable  conductors  has  not  been  fully 
explored  in  4D  printing  and  soft  machines.  This  work  studies  the  effect  of  thermal  cure  conditions  on 
conductivity  and  electro-mechanical  behaviors  of  silver  ink  by  the  direct  ink  write  (DIW)  printing  approach. 

This  work  is  published  in  "Mu,  Q,  Dunn,  C.K.  Wang,  L.,  Dunn,  M.L.,  Qi,  H.J.,  Wang,  T.,  2016.  Thermal 
cure  effects  on  electromechanical  properties  of  conductive  wires  by  direct  ink  write  for  4D  printing  and 
soft  machines,  submitted." 

Our  studies  showed  that  that  the  electro-mechanical  properties  of  silver  wires  can  be  tailored  by 
controlling  cure  time  and  cure  temperature  to  achieve  conductivity  as  well  as  stretchability.  For  the  silver 
NP  ink  we  used  in  the  experiments,  silver  wires  cured  at  80°C  for  10-30min  have  conductivity  >1%  bulk 
silver,  Young's  modulus  <100MPa,  yield  strain  ~9%,  and  can  retain  conductivity  up  to  300%  strain.  In 
addition,  under  stress  controlled  cyclic  loading/unloading  conditions,  the  resistance  of  these  wires  is  only 
about  1.3  times  the  initial  value  after  the  100th  repeat  cycle  (7.6%  maximum  strain  in  the  first  cycle).  Silver 
wires  cured  at  120°C  for  10-20min  are  more  sensitive  to  strain  and  have  a  yield  strain  of  around  6%.  These 
properties  indicate  that  the  silver  ink  can  be  used  to  fabricate  stretchable  electrodes  and  flex  sensors. 
Using  the  DIW  fabrication  method,  we  printed  silver  ink  patterns  on  the  surface  of  3D  printed  polymer 
parts,  with  the  future  goal  of  constructing  fully  3D  printed  arbitrarily  formed  soft  and  stretchable  devices 
and  of  applying  them  to  4D  printing.  We  demonstrated  a  fully  printed  functional  soft-matter  sensor  and 
a  circuit  element  that  can  be  stretched. 

Materials  and  Method:  The  silver  nanoparticle  ink  used  in  this  study  was  DuPont  ME603,  which 
was  donated  by  DuPont  Inc.  and  has  a  solid  content  of  49%-53%,  viscosity  of  15-35mP*s  and  sheet 
resistivity  <200mQ/sq/mil.  The  manufacturer  recommends  the  curing  condition  of  120°C  -140°C  in  a  static 
box  oven  for  lOmin.  In  this  work,  we  DIW  printed  the  ink  onto  a  3D  printed  elastomer  substrate,  namely 
TangoBlack,  one  of  the  model  materials  available  in  the  material  library  of  the  polyjet  3D  printer  Objet 
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Connex  260  (Stratasys  Inc.,  Edina,  MN,  USA).  The  reason  for  choosing  a  3D  printed  substrate  is  that  our 
work  is  intended  to  combine  3D  printed  materials  with  DIW  of  functional  inks,  and  TangoBlack  is  a  material 
used  in  our  previous  work  on  4D  printing.  In  addition,  we  found  that  when  the  silver  ink  was  printed  onto 
a  soft  substrate  such  as  an  elastomer,  it  could  be  easily  peeled  off  for  mechanical  testing. 

Direct  Ink  Write  setup:  We  printed  the  silver  ink  onto  a  3D  printed  elastomer  substrate  using  a 
DIW  experimental  setup,  which  is  shown  in  Figure  2.6.1-la.  The  printing  platform  moves  the  substrate 
along  the  y-direction,  while  a  pressure  controlled  syringe  moves  in  the  x-  and  z-  directions.  The  pressure 
in  the  syringe  is  controlled  by  a  pressure  regulator  (Ultimus  V,  Nordson  EFD,  East  Providence,  Rl,  USA). 
The  syringe  was  attached  to  a  nozzle  with  0.61mm  inner  diameter.  With  a  pressure  of  15psi  and  a  printing 
speed  of  3mm/s,  we  were  able  to  print  the  silver  ink  with  consistent  cross-sectional  areas. 

Thermal  curing:  Once  the  silver  ink  was  printed,  it  was  cured  in  a  convection  oven  (Model  13-247- 
750G,  Fisher  Scientific,  Pittsburgh,  PA,  USA)  at  different  temperatures  and  curing  times.  To  record  the 
resistance,  change  during  curing,  we  fabricated  two  copper  electrodes  that  connected  the  printed  lines 
to  a  multimeter  (Keithley2100,  Keithley  Instruments,  Cleveland,  OH,  USA),  which  was  controlled  by  a 
LabVIEW  data  acquisition  program  and  sampled  at  0.5-5Hz. 

Measurement  of  cross-sectional  area:  To  calculate  the  resistivity  and  modulus  of  printed  silver 
wires  cured  under  different  conditions,  the  cross  sectional  area  of  these  wires  should  be  known.  Cross- 
sectional  images  of  the  cured  silver  wires  were  taken  with  a  trinocular  microscope  (Model  89404-886, 
VWR,  Suwanee,  GA,  USA).  2.6.1-lc  shows  a  typical  cross-section  of  the  cured  ink.  The  cross-sectional  area 
was  then  calculated  by  using  the  open  source  image  processing  software  Image! 

Stress-strain  and  cyclic  loading/unloading  tests:  After  curing  and  measuring  the  cross-sectional 
area,  the  silver  wires  were  carefully  peeled  off  from  the  substrates  using  a  pair  of  tweezers.  The  electrical 
and  mechanical  properties  were  then  studied.  We  used  a  dynamic  mechanical  analysis  (DMA)  tester 
(Model  Q800,  TA  Instruments,  New  Castle,  DE,  USA)  and  a  multimeter  to  measure  the  electro-mechanical 
properties  of  cured  freestanding  wires.  For  the  stress-strain  tests,  the  strain  control  mode  was  used  with 
a  strain  rate  of  10%  /min.  For  the  cyclic  loaded/unloaded  tests,  the  stress  control  mode  was  used  with  a 
stretch/release  speed  of  25MPa/min  (corresponding  to  a  strain  rate  of  0.5%/s  -2%/s),  and  a  holding  time 
of  0.5min  between  cycles.  The  resistance  change  during  cyclic  testing  was  also  recorded.  All 
electromechanical  property  tests  were  performed  at  room  temperature  (22°C). 


Figure  2. 6. 1-1.  (a)  A  schematic  graph  of  the  DIW  setup  where  the  syringe  can  move  in  the  x-  and  z- 
directions  and  the  motion  stage  can  move  in  the  y-direction.  A  pressure  regulator  is  used  to  control  the 
pressure  applied  to  the  syringe,  (b)  A  photo  of  the  DIW  experimental  setup,  (c)  A  microscope  image  of  a 
cross  section  of  written  silver  wire  on  the  elastomer  substrate. 
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Resistivity:  For  all  resistance  measurements,  a  two-point  probe  method  was  used  to  calculate  the  bulk 
resistivity  and  its  change  throughout  curing. 


Results  and  discuss  -  Effect  of  thermal  cure  conditions:  In  order  to  choose  proper  curing  temperatures, 
both  TGA  and  DSC  tests  were  conducted  to  study  the  mass  change  during  curing.  The  TGA  and  DSC  curves 
are  shown  in  Figure  2.6. 1-2.  It  can  be  seen  that  the  mass  of  the  sample  decreases  as  the  temperature 
increases,  settling  at  ~54%  of  the  initial  mass.  This  value  confirms  the  removal  of  solvents  and  additives 
and  agrees  with  the  solid  content  (49%-53%)  of  the  ink  according  to  the  manufacturer.  The  mass  reduction 
can  be  divided  into  three  stages,  which  are  shown  by  the  slope  change  in  the  TGA  curve,  and  which  can 
be  identified  from  the  three  peaks  in  the  derivative  thermolgravimetry  (DTG)  curve.  In  the  first  stage,  the 
sample's  mass  quickly  drops  26%  between  47°C  and  177°C,  with  an  endothermic  reaction  observed  from 
the  DSC  curve,  which  may  be  caused  by  the  evaporation  of  solvents.  In  the  second  stage,  the  mass  sharply 
drops  another  11%  near  ~200°C,  where  the  peak  on  the  DSC  test  indicates  decomposition  of  organic 
additives  or  removal  of  residual  solvents.  Finally,  in  the  third  stage,  at  around  440°C  there  is  another  mass 
drop,  which  may  be  caused  by  decomposition  of  other  organic/inorganic  additives. 

In  situ  resistance  change  during  curing:  The  TGA/DSC  curves  show  that  as  the  temperature 
increases,  the  solvent  in  the  ink  evaporates  at  different  temperatures.  Thus,  we  cured  the  ink  at  different 
temperatures  in  the  range  of  80°C  -180°C  for  a  relatively  long  time  (60min)  to  detect  resistance  changes. 
We  compared  the  curing  processes  at  four  different  temperatures:  80°C,  100°C,  120°C,  and  140°C  (Figure 
2.6.1-3a).  To  reduce  the  variations  between  individual  samples,  we  normalized  the  resistance  during 
curing  by  using  individual  samples'  respective  Reomin,  the  resistance  measured  at  60min.  The  infinitely  large 
resistance  at  the  onset  of  curing  is  because  the  silver  NPs  are  surrounded  and  separated  by  solvent  and 
organic  shell,  inhibiting  the  formation  of  a  conductive  path.  The  resistance  changes  show  that  the  solvent 
evaporates  gradually  as  a  function  of  curing  time,  where  longer  curing  times  contribute  to  the  drop  in 
resistivity  and  the  formation  of  conductive  paths  within  wires.  For  the  120°C  and  140°C  cured  samples, 
the  resistance  is  nearly  saturated  within  the  first  lOmin,  which  agrees  with  manufacturer  suggested  curing 
conditions,  while  the  80°C  and  100°C  cured  samples  show  slower  curing  rates,  but  the  resistivity  finally 
becomes  constant  after  curing  for  40min-60min. 


0  100  200  300  400  500  600 
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Figure  2. 6. 1-2.  TGA  (red  line)  and  DSC  (black  line)  measurements  of  the  ink  during  curing:  Heated  from 
40.00°C  to  600.00°C  at  10.00°C  /min.  The  specimen  mass:  18.343  mg  for  TGA,  and  13.790  mg  for  DSC.  The 
derivative  thermolgravimetry  (DTG)  curve  is  also  shown  in  black  color. 
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Figure  2. 6. 1-3.  (a)  Resistance  change  during  curing  at  80°C,  100°C,  120°C,  140°C  for  about  60min.  (b) 
Calculated  resistivity  cured  at  different  temperature  for  60min  and  compared  with  bulk  silver. 

To  calculate  the  resistivity,  optical  microscope  images  of  the  wire's  cross  section  were  taken  to 
measure  the  cross  sectional  area  (see  Figure  2.6.1-lc  for  a  typical  image).  Under  the  aforementioned 
operation  conditions  (15psi  pressure,  3mm/s  speed),  the  wires  had  an  average  width  of  1.07mm,  and  a 
cross-sectional  area  of  0.062mm2.  Using  this  measurement,  the  resistivity  was  calculated  and  is  shown  in 
Figure  2.6.1-3b  for  curing  60min  at  different  temperatures.  The  resistivity  of  bulk  silver  is  also  shown  for 
comparison.  With  curing  temperatures  greater  than  140°C,  the  effect  of  temperature  on  resistivity  is  not 
obvious,  and  the  saturated  resistivity  is  ~8  times  higher  than  that  of  bulk  silver.  In  the  lower  curing 
temperature  range  (80°C  -120°C),  the  resistivity  drops  significantly  as  the  temperature  increases. 
Interestingly,  even  at  80°C,  the  ink  can  be  conductive  after  being  cured  for  60min,  and  the  resistivity  is 
only  5-6  times  higher  than  those  cured  at  140°C.  In  addition,  Figure  2.6.1-3a  shows  that  the  resistance 
becomes  almost  constant  after  40min-60min  (even  at  80°C),  indicating  a  wire  with  stable  conductivity  is 
formed. 

Microstructure  characterization:  SEM  images  (Figure  2. 6. 1-4)  show  the  morphology  and  particle  packing 
for  the  ink  cured  under  different  conditions.  For  comparison,  images  from  the  air-dried  ink  are  also 
presented.  The  overall  view  (image  in  the  middle)  of  the  air-dried  ink  shows  separated  small  islands  with 
different  shapes  and  diameters  without  a  clear  interconnected  path.  A  zoomed-in  image  (the  left  image) 
of  the  connected  area  shows  that  the  silver  nanoparticles  of  ~lpm  diameter  are  nearly  separated  from 
each  other;  for  the  unpacked  area,  small  silver  particles  having  diameters  of  ~100nm  are  separated  from 
each  other  without  forming  a  conductive  path.  In  samples  cured  at  80°C  for  30min,  we  see  the  formation 
of  conductive  paths.  As  shown  in  the  image  in  the  middle,  connected  islands  are  forming,  with  some 
amount  of  holes.  The  zoom-in  of  the  connected  area  (left  image)  shows  fully  packed  silver  particles  and 
inter-particle  boundaries  can  be  seen  but  not  very  clearly.  For  the  unpacked  area  (right  image),  we  see 
large  particles  that  are  close  to  each  other  (as  compared  to  the  air-dried  samples).  When  cured  at  180°C 
for  30min,  some  pores  exist,  which  might  be  responsible  for  the  relatively  large  resistivity  as  compared  to 
bulk  silver,  but  large  connected  regions  are  also  obvious.  From  the  zoom-in  of  the  connected  region  (left 
image),  we  can  see  better-connected  particles;  even  in  the  not-fully  connected-region  (right  image),  silver 
NPs  are  physically  in  contact. 

Electro-mechanical  properties!  We  conducted  tensile  tests  on  freestanding  silver  wires.  Three  samples 
were  tested  for  each  curing  condition.  Figure  2. 6. 1-5  shows  results  from  two  representative  samples 
subject  to  two  different  curing  conditions  (at  80°C  for  20min  and  at  140°C  for  20min,  respectively).  The 
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sample  cured  at  a  higher  temperature  exhibits  a  much  higher  yield  stress  and  lower  yield  strain  than  does 
the  lower  temperature  cured  sample.  The  resistance  of  the  high  temperature  cured  sample  is  more 
sensitive  to  strain  during  stretching  than  that  of  the  low  temperature  cured  one.  In  addition,  surprisingly, 
the  silver  wire  can  be  stretched  significantly.  For  the  sample  cured  at  140°C  for  20min,  the  cured  wire  can 
be  stretched  ~150%,  but  with  significantly  increased  resistance;  the  sample  cured  at  80°C  for  20min  cured 
can  be  stretched  300%  with  an  increase  of  resistance  of  ~  400  times. 


Figure  2. 6.1-4.  Comparison  of  SEM  image  of  surface  morphology:  as  dried,  cured  for  30  min  at  80°C  and 
cured  for  30  min  at  180°C. 


Strain(%) 

Figure  2. 6. 1-5.  Stress-strain  curves  of  printed  silver  wires  and  normalized  resistance  change  when 
stretched,  taking  80°C  20min  and  140°C  20min  cured  samples  as  examples.  The  inserted  image  shows  how 
to  obtain  yield  strain  and  stress  from  the  stress-strain  curve  of  the  140°C  20min  cured  sample. 
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For  applications  in  4D  printing  and  soft  machines,  the  conductive  materials  should  have  1)  high 
electrical  conductivity  to  achieve  high-efficiency  for  the  devices,  2)  compliance  that  matches  with  the 
substrate,  3)  high  elasticity,  and  4)  an  almost  constant  electrical  conductivity  under  applied  strain.  In  the 
following,  we  focus  on  the  conductivity,  modulus,  yield  strain,  and  elongation  of  the  DIW  printed  silver 
wires. 

Conductivity:  Based  on  the  curing  condition  study,  we  chose  to  cure  the  silver  NPs  in  the  oven  at  six 
different  temperatures:  80°C,  100°C,  120°C,  140°C,  160°C  and  180°C,  with  curing  durations  of  lOmin, 
20min,  30min,  40min,  50min  and  60min;  the  results  are  summarized  in  Figure  2.6.1-6a.  In  general,  the 
resistivity  decreases  as  the  curing  temperature  and  curing  time  increase.  For  the  range  of  these  curing 

conditions,  the  resistivity  is  between  150.9pQcm  -10.2pQcm,  which  corresponds  to  a  resistivity  of  94 
and  6.4  times  that  of  bulk  silver  (1.6  pQcm  [14]),  respectively.  Since  the  3D  printed  materials  from  the  3D 
printer  can  sustain  a  high  temperature  of  ~120°C,  we  chose  the  silver  wire  that  was  cured  under  the  mild 
temperature  (lower  than  120°C)  to  serve  as  conducting  wires  for  4D  printing  and  soft  machine 
applications. 

Young's  moduli:  Figure  2.6. 1-7  shows  the  evolution  of  Young's  moduli  as  a  function  of  curing  temperature 
and  curing  time.  Here,  the  Young's  modulus  is  determined  from  the  initial  slope  of  the  stress-strain  curve 
(~0.5%).  It  is  clear  that  a  higher  temperature,  or  a  longer  curing  time,  leads  to  a  higher  Young's  modulus. 
This  is  consistent  with  the  SEM  observations  that  higher  curing  temperature  yields  to  more  compacted 
and  better  connected  silver  NPs.  The  Young's  moduli,  however,  are  lower  than  that  of  bulk  silver.  From 
Figure  2.6. 1-7,  the  highest  modulus  is  ~800MPa,  under  the  curing  condition  of  140°C  for  60min.  For  bulk 
silver,  the  modulus  of  electron-beam  deposited  free  stranding  Ag  films  is  roughly  63GPa,  and  the 
calculated  E  <111>  is  82.2  GPa  [15];  inkjet  printed  silver  film  after  being  sintered  at  250°C  has  a  modulus 
around  50GPa  [16].  It  is  also  interesting  to  note  that  for  the  80°C  curing  temperature  we  observe  no 
significant  change  in  modulus  as  a  function  of  curing  time,  even  though  this  temperature  yields  a 
conductive  line  after  curing.  In  our  experiments,  we  found  that  curing  temperatures  higher  than  140°C 
resulted  in  wires  that  were  very  stiff,  which  significantly  limited  their  ability  to  stretch.  Therefore,  in  the 
following,  samples  with  curing  temperatures  below  140°C  were  used  for  the  electromechanical  property 
studies. 


T(°C) 


181.40  143.52  111.47  83.99 


(a)  (b) 

Figure  2.6. 1-6.  (a)  Conductivity  variation  of  the  silver  conductive  wires  as  a  function  of  cure  temperature 
and  cure  time,  (b)  Arrhenius  plot  of  electrical  conductivity. 
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Figure  2.6. 1-7.  Modulus  as  function  of  cure  time  and  cure  temperature. 


(a)  (b) 

Figure  2.6. 1-8.  Yielding  (point  of  intersection  of  tangents)  behavior  as  a  function  of  cure  time  and 
temperature  (a)  yielding  strain;  (b)  yielding  stress. 

Yield  strain:  We  determine  the  yielding  point  from  the  intersection  of  tangents  of  the  stress-strain  curves 
to  obtain  the  yielding  stain  and  stress.  The  yielding  strains  and  yielding  stresses  are  summarized  in  Figure 

2.6. 1- 8.  The  yielding  strain  decreases  with  an  increase  in  curing  time  and  temperature,  while  the  yielding 
stress  increases  with  an  increase  in  curing  time  and  temperature.  In  general,  the  cured  silver  wires  have 
yielding  strains  between  2.0%  and  9.4%,  which  is  larger  than  freestanding  metal  film  whose  stretch  is  only 
1%-2%[17]. 

Elongation:  As  in  Figure  2.6. 1-5  shows,  the  cured  silver  wires  can  withstand  a  large  elongation.  During  the 
experiments,  the  resistance  increased  into  the  MQ  range  when  the  wire  was  subjected  to  large  strains 
(>300%).  Here,  for  the  sake  of  analysis,  we  define  the  ultimate  elongation  as  the  strain  when  the 
normalized  resistance  is  increased  by  100  times  (R/Ro  <100)  or  the  silver  wire  breaks  and  conductivity  is 
lost,  whichever  occurs  first.  Figure  2.6. 1-9  shows  the  ultimate  elongation  of  wires  cured  under  different 
conditions.  In  general,  we  can  obtain  ~100%  ultimate  elongation  for  most  curing  conditions,  which  is  much 
larger  than  sputtered  metal  films,  which  have  a  maximum  elongation  of  around  50%  on  elastic 
substrate[18,  19].  This  value  is  also  larger  than  screen  printed  silver  ink,  which  is  about  50%[20].  Figure 

2.6. 1- 9  also  shows  lower  temperature  cured  silver  wire  has  a  longer  elongation.  In  addition,  20min  curing 
time  gives  the  largest  ultimate  elongation  for  all  the  curing  temperatures. 
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Cure  time  (min) 

Figure  2.6. 1-9.  Ultimate  elongation  (R/R0  <100  or  sample  crack)  as  a  function  of  cure  temperature  and 
cure  time. 


Figure  2.6.1-10  shows  that  the  silver  ink  in  this  study  is  capable  of  large  elongation.  Therefore,  if 
it  is  deposited  on  an  elastic  substrate  it  may  achieve  high  elastic  deformation  by  embedded  printing[21], 
encapsulation[22],  and  other  methods  that  have  potential  for  stretchable  electrodes  application.  In 
addition,  we  found  that  for  the  silver  wires  we  printed,  the  cure  temperature  of  80°C  -120°C  for  10-30min 
is  suitable  for  4D  printing  and  soft  machine  usage,  which  typically  require  moderate  to  large  deformation 
of  the  substrate  materials 

Examples  of  potential  applications:  Based  on  the  thermal  cure  condition  as  well  as  the 
electromechanical  properties  study,  we  found  that  as  the  curing  temperature  and  curing  time  increased 
the  cured  wires  became  more  conductive  but  less  stretchable.  The  silver  wires  cured  at  80°C  for  10-30min 
are  stretchable,  and  thus  can  be  used  for  stretchable  applications;  the  silver  wires  cured  at  120°C  for  10- 
20min  have  good  conductivity  and  are  sensitive  to  applied  strain;  they  are  thus  good  for  sensing.  In  the 
following,  we  demonstrate  two  potential  applications. 

The  printed  ink  can  be  embedded  into  polyJet  3D  printed  parts,  resulting  in  the  fabrication  of  fully 
3D  printed,  soft  and  stretchable  devices.  To  fabricate  such  devices,  we  used  a  sequential  process,  where 
we  printed  a  substrate  using  the  3D  printer,  followed  by  DIW  printing,  curing  the  ink  and  then 
encapsulating  the  ink  by  subsequent  3D  printed  layers.  This  process  can  be  repeated  to  create  multiple 
layers  of  embedded  ink.  As  a  demonstration,  a  stretchable  LED  circuit  was  printed,  shown  in  Figure  2.6.1- 

10,  where  eight  ink  lines  were  patterned  inside  a  60mm  long  strip  of  TangoBlack  connecting  four  LEDs  in 
parallel,  mimicking  a  typical  ribbon  cable  with  multiple  conducting  wires.  The  ink  wires  were  cured  at  80°C 
for  20  minutes  in  order  to  yield  highly  stretchable  lines  while  still  maintaining  a  high  conductivity  to 
illuminate  the  LEDs  at  a  reasonable  power  output  (~500mW).  It  was  possible  to  stretch  the  cable  by  as 
much  as  45%  during  testing.  The  online  video  shows  the  stretching  process. 

A  flexible  sensor  was  also  printed.  Here,  a  single  layer  serpentine  strain  gauge  was  printed  inside 
of  the  elastomer  ring.  The  ink  line  was  cured  at  120°C  for  20  minutes,  a  more  suitable  curing  condition  for 
sensor  applications.  The  fabricated  sensor  and  the  typical  operating  response  can  be  seen  in  Figure  2.6.1- 

11.  The  sensor  was  tested  in  two  bending  configurations,  1)  straight  and  2)  partially  bent.  When  the  finger 
bends,  it  bends  the  top  face  of  the  ring  and  stretches  the  conductive  wire,  causing  an  increase  of  the 
resistance.  When  the  finger  returns  to  the  straight  position,  the  resistance  drops.  The  video  in  the  online 
material  shows  the  real  time  resistance  change  as  the  finger  bends  and  returns. 
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Figure  2.6.1-10.  Fully  printed  stretchable  cable  by  embedded  silver  wire  into  polyJet  3D  printed  parts:  (a) 
The  design  of  the  cable,  (b)  The  cross-section  of  the  stretchable  part  of  the  cable,  (c)  The  images  of  the 
printed  part  during  the  printing  process,  (d)  Embedded  silver  wires  in  elastomer  as  circuits  during 
stretching.  Images  of  the  printed  cable  under  tensile  strain  of  (  1)0%  and  (II)  40%  and  (HI)  under  folding. 


Figure  2.6.1-11.  Detection  of  finger  motion  with  all  printed  stretchable  sensor,  (a)  Schematic  and  printed 
flex  sensor,  (b)  Images  of  a  bending  sensor  on  a  finger-ring  mounted  on  the  forefinger,  in  straight  (I)  and 
bend  (II)  positions,  (c)  Change  in  resistance  of  a  representative  gauge  during  four  bending  cycles. 
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2.7  Design  Optimization 

To  develop,  implement,  and  study  our  topology  optimization  framework,  we  first  identified  and 
implemented  an  accurate  and  efficient  solid  shell  finite  element  formulation  in  our  in-house  design 
optimization  software  platform.  The  solid  shell  formulation  allows  modeling  and  predicting  the 
mechanical  phenomena  in  composite  plates  and  shells  with  complex  material  arrangements.  We  further 
developed  and  studied  a  density  method  for  optimizing  active  composites.  This  study  revealed  significant 
shortcomings  of  the  density  approach  for  the  design  of  PACS.  Therefore,  we  focused  on  developing  and 
applying  a  LSM-XFEM  approach  for  finding  the  optimal  arrangement  of  active  materials  within  an  inactive 
phase.  To  facilitate  optimizing  PACS,  we  develop  first  a  simplified  PAC  model  amendable  to  iterative 
design  optimization  procedures.  We  refined  and  applied  our  LSM-XFEM  approach  for  finding  the  optimal 
arrangement  of  active  materials  within  an  inactive  phase.  These  studies  involved  printing  samples  of 
optimized  PACs  and  verifying  their  behavior  experimentally.  While  the  simplified  PAC  model  captured  the 
basic  shape  memory  polymer  phenomena,  it  relied  on  a  small  strain  assumption,  which  severely  limits  the 
applicability  of  the  optimization  approach.  Therefore,  we  extended  our  simplified  PAC  model  to  account 
for  finite  strains.  We  further  refined  our  models  and  LSM-XFEM  optimization  approach  for  finding  the 
optimal  arrangement  of  heterogeneous  materials  with  imperfect  bonding  between  the  material  phases. 
A  foundational  study  on  design  sensitivities  with  the  Heaviside  enrichment  strategy  was  a  significant 
byproduct  of  this  work. 

2.7.1  Solid  Shell  Finite  Element  Formulation  and  Implementation 

Conventional  finite  element  formulations  of  plate  and  shell  models  are  based  on  a  degeneration  principle, 
i.e.  the  behavior  in  thickness  direction  is  modeled  explicitly  and  the  structural  response  is  represented  by 
the  deformations  of  the  mid-plane.  While  computationally  efficient  and  sufficiently  accurate  for  regular 
material  arrangements,  this  approach  is  inadequate  for  complex  material  layouts  and  complicates 
modeling  complex  3D  constitutive  models.  These  issues  can  be  mitigated  by  using  solid  shell  elements 
which  are  built  upon  a  3D  formulation  of  the  mechanical  model.  However,  standard  solid  elements  suffer 
from  various  locking  phenomena  and  require  highly  resolved  meshes,  in  particular  in  thickness  direction, 
to  predict  accurately  the  response  of  thin  structures.  Solid  shell  elements  mitigate  this  issue  by  enriching 
the  strain  approximation  via  Assumed  Natural  Strain  (ANS),  Enhanced  Assumed  Strain  (EAS)  concepts,  and 
tuned  numerical  integration  procedures. 


Figure  2.7. 1-1.  Benchmark  example  for  solid  shell  element  implementation:  Cylinder  buckling  under 
tensile  forces. 
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displacements  at  points  A,B,C 


Figure  2. 7. 1-2:  Load-displacement  curves  for  points  A,  B,  and  C  on  cylinder;  see  Figure  2.7. 1-1. 

To  enable  the  modeling  of  complex  material  layouts  by  finite  elements,  we  identified  and  implemented  a 
solid  shell  element  formulation  that  does  not  rely  on  simplifying  assumption  of  the  material  layout  and 
that  can  be  integrated  into  an  XFEM  formulation.  After  studying  the  literature  and  implementing  6  of  the 
most  promising  formulations  into  a  MATLAB  prototype  code,  we  identified  the  formulation  of  Schwarze 
and  Reese  [23,  24]  as  best  suited  for  meeting  the  needs  of  this  project.  This  solid  shell  formulation  was 
implemented  in  our  in-house  finite  element  &  optimization  software  platform,  verified  against  a 
comprehensive  suite  of  benchmark  problems,  and  has  been  used  in  the  design  optimization  studies 
described  below.  Figure  2.7. 1-1  shows  one  of  the  benchmark  problems  studied.  A  cylinder  is  subject  to 
two  vertical  forces  acting  in  opposite  directions.  The  cylinder  undergoes  a  highly  nonlinear  behavior  as 
the  load-displacement  curves  for  the  points  A,  B,  and  C  illustrate;  see  Figure  2. 7. 1-2.  Our  implementation 
matches  perfectly  the  results  reported  in  [24]. 


2.7.2  Optimization  of  Printed  Active  Composites  by  the  Density  Approach 

We  first  explored  the  capabilities  to  optimize  the  material  layout  of  printed  active  composites  with  the 
density  approach.  We  describe  the  response  of  the  composite  via  a  homogenized  macroscopic  material 
model  assuming  small  strains.  The  material  properties  are  derived  from  a  microscopic  model  using  Mori- 
Tanka's  homogenization  approach  [25,  26].  The  micro-scale  material  architecture  consists  of 
unidirectional  SMP  fibers  embedded  in  an  isotropic  inactive  matrix.  The  shape  memory  effect  of  the  active 
fibers  is  modeled  via  incompressible  eigenstrains  which  approximate  residual  strains  due  to  the  training 
load  via  the  magnitude  and  orientation  of  the  eigenstrains.  The  fiber  volume  fraction  and  the  fiber 
orientations  were  considered  optimization  variables  and  their  optimum  spatial  distribution  determined  in 
the  optimization  process. 

In  our  design  studies,  we  determined  the  optimal  design  of  printed  active  composites  such  that 
in  the  actuated  stage,  i.e.  after  releasing  the  training  loads  at  TL,  the  deformed  shape  of  a  composite  plate 
matches  a  target  shape.  An  example  is  shown  in  Figure  2. 7. 2-1.  The  objective  of  the  problem  is  to  find  the 
optimum  distributions  of  fiber  volume  fraction  and  in-plane  fiber  orientation,  i.e.  in  the  x-y  plane,  such 
that  a  square  composite  plate  assumes  a  parabolic  shape  upon  actuation,  i.e.  due  to  the  shape  memory 
behavior  of  the  fibers.  The  design  domain  is  restricted  to  the  lower  half  of  the  plate.  A  residual  fiber  strain 
in  the  x-direction  of  4%  and  in  the  y-  and  z-direction  of -2%  is  assumed  to  capture  the  training  loads.  Figure 
2. 7. 2-1  shows  the  optimum  distributions  of  fiber  volume  fraction  and  in-plane  fiber  orientation,  as  well  as 
the  deformed  and  target  shapes.  Note  the  fiber  orientations  are  visualized  by  drawing  "streamlines"  along 
the  fiber  directions.  The  optimized  composite  layout  results  in  an  actuated  shape  that  matches  well  the 
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Fiber  volume  fraction 


Fiber  volume  fraction 
with  "stream  lines"  along  fiber  direction 


Figure  2. 7. 2-1.  Optimization  of  fiber  volume  fraction  and  fiber  orientation  for  shape  matching. 


target  shape.  While  a  constraint  on  the  maximum  total  fiber  volume  fraction  of  50%  was  imposed,  this 
constraint  was  not  active  in  the  final  design.  Studies  on  different  target  shapes  were  performed  and  also 
showed  that  a  broad  range  of  target  shapes  can  be  matched  with  great  accuracy  through  optimizing  the 
fiber  volume  fraction  and  fiber  orientation. 

We  studied  methods  for  converting  the  optimized  fiber  volume  fraction  and  fiber  orientations 
into  a  printable  design.  To  this  end,  we  represented  individual  fibers  by  tubes  which  were  generated  as 
streamlines  along  the  principal  stiffness  directions.  The  tube  radius  was  defined  as  a  function  of  the  local 
fiber  volume  fraction.  The  parameters  of  this  function  and  the  spacing  of  the  fibers  were  tuned  manually 
to  match  the  total  fiber  volume  fraction  in  the  plate.  We  then  generated  surface  meshes  for  the 
tube/matrix  interfaces  which  were  used  by  the  printing  software  of  the  Object  260  printer  to  fabricate 
the  composite;  see  Figure  2. 7. 2-2.  However,  preliminary  experimental  studies  suggest  that  this  heuristic 
approach  is  prone  to  introducing  significant  discrepancies  between  the  mechanical  model  used  for 
optimization  and  the  model  extracted  for  printing. 


meshed  fibers  (STL  format) 


printed  design  (fibers  -  black) 


Figure  2. 7. 2-2.  Realization  of  optimized  active  composite. 
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2.7.3  Simplified  PAC  Model  for  Design  Optimization 

To  accurately  model  the  response  of  PACs,  in  general  large  deformations  and  a  complex  thermo¬ 
mechanical  constitutive  behavior  need  to  be  considered,  leading  to  a  nonlinear  non-conservative  multi¬ 
physics  model.  Constitutive  models  for  accurately  describing  the  SM  cycle  of  amorphous  shape  memory 
polymers  (SMPs)  have  been  developed  in  the  first  year  of  this  project  and  are  described  in  the  previous 
project  report.  Considering  these  models  within  the  topology  optimization  process  significantly  increases 
the  algorithmic  complexity  and  the  computational  cost.  To  create  an  analysis  approach  more  amenable 
for  design  in  terms  of  computational  speed  and  ease  of  implementation  but  with  acceptable  fidelity,  we 
developed  two  simplified  mechanical  models.  We  studied  the  accuracy  of  these  PAC  models  against 
experimental  results,  and  compared  the  optimization  results  obtained  using  these  models.  This  work  was 
published  in  K.  Maute,  A.  Tkachuk,  J.  Wu,  H.  J.  Qi,  Z.  Ding ,  and  M.L  Dunn.  2015.  Level  Set  Topology 
Optimization  of  Printed  Active  Composites.  ASME  Journal  of  Mechanical  Design ,  137: 111704-1,  2015. 

Both  simplified  PAC  models  account  for  stress  relaxation  and  the  resulting  residual  stresses  via 
eigenstrains  that  can  vary  with  position.  These  models  lump  the  complex  nonlinear  thermo-mechanical 
phenomena  occurring  during  a  thermo-mechanical  cycle  into  an  effective  inelastic  strain,  i.e.  eigenstrain. 
Both  PAC  models  assume  infinitesimal  strains  and  a  linear  material  behavior,  but  they  differ  in  the 
approximation  of  the  eigenstrains  in  the  active  phase.  The  more  accurate  model  predicts  the  eigenstrains 
by  resolving  the  strain  field  in  the  PAC  at  the  training  stage  due  to  the  programming  loads  at  temperatures 
above  the  glass  transition  temperature.  The  simpler  model  assumes  a  uniform  eigenstrain  in  the  active 
phase  which  is  predicted  analytically.  Subsequently,  we  first  describe  the  more  accurate  model  which  is 
then  simplified  to  the  model  with  uniform  eigenstrains.  Both  models  were  compared  against  experimental 
data. 

We  consider  the  response  of  the  PAC  at  the  training  stage  subject  to  programming  loads,  referred 
to  stage  (i),  and  after  the  programming  loads  are  removed,  referred  to  stage  (ii).  The  effects  of  the  phase 
transition  of  the  active  phase  are  modeled  via  eigenstrains,  which  are  computed  in  stage  (i)  and  applied 
in  stage  (ii),  and  a  change  in  Young's  modulus.  The  residuals  of  the  static  equilibrium  equations  in  both 
stages  are  written  in  weak  form  as  follows: 

%  =  t(L^  dQ-j^Su'F,'  rfr)  =  0,  (2.73-1) 

t  dn-jt;SufF,'  rfr)=o,  (2.73-2) 

with 

up  =UP  onf£,  up  =  Upo nf£,  (2.7.3-3) 

where  (“)  denotes  quantities  in  stage  (i)  and  (")  in  stage  (ii).  Passive  and  active  phases  are  denoted 
with  superscripts  "1"  and  "2".  The  displacement  vector,  the  strain  tensor,  and  the  stress  tensor  are 
denoted  by  up ,  sp  and  op  respectively.  External  forces,  Fp ,  act  along  the  boundary  Yp  and 

prescribed  displacements,  Up  ,  are  applied  at  the  boundary]"^  .  Note  the  boundary  conditions  typically 

differ  between  stages  (i)  and  (ii).  The  material  behavior  at  the  two  stages  is  described  by  the  following 
constitutive  equations: 


(2. 7.3-4) 
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(2.73-5) 


<  = 


where  CJjkl ,  fp,  and  ap  denote  the  elastic  tensor,  the  strain  fixity  and  the  coefficient  of  thermal 

expansion,  respectively.  The  Young's  modulus  of  the  active  phase  increases  as  the  material  transitions 
from  a  rubbery  state  in  stage  (i)  to  a  glassy  state  in  stage  (ii).  The  elastic  tensor  of  the  passive  phase  is 
assumed  constant.  The  strains  of  the  programming  state  (i)  are  introduced  as  eigenstrains  in  the 

constitutive  equation  for  stage  (ii).  For  the  passive  phase  the  fixity  is  assumed  negligible,  i.e.  fx  —  0.0  .  In 

the  following,  we  ignore  thermal  strains  as  the  mismatch  in  ap  is  insignificant  for  the  materials 
considered  in  this  study. 

Subsequently,  we  refer  to  the  model  described  above  as  two-stage  model.  This  model  can  be 
further  simplified  by  approximating  the  strains  in  the  training  stage  (i)  by  a  uniform  strain  distribution. 

Assuming  a  homogeneous  material  and  for  simple  geometries  and  programming  loads,  ~ep  can  be 

approximated  analytically  via  kinematic  compatibility  considerations.  As  before,  the  fixity  in  the  passive 
phase  is  assumed  to  vanish.  Note  that  this  simplification  requires  only  modeling  the  static  equilibrium  of 
stage  (ii).  We  refer  to  this  model  as  "single  stage"  model. 

The  PAC  models  introduced  above  are  built  upon  several  simplifying  assumptions.  To  investigate 
the  ability  of  these  models  to  capture  the  fundamental  characteristics  of  the  PAC  response  we  consider 
the  simple  configuration  of  Figure  2.7. 3-1,  which  was  studied  experimentally  [29].  The  active  fibers  are 
made  of  Shore95  and  exhibit  a  SM  behavior;  the  passive  phase  is  Tango  black+.  In  the  programming  stage 
(i),  the  composite  plate  is  subject  to  prescribed  displacements  ±0.56^  at  the  front  and  rear  faces.  In  the 


Level  set  field:  load  direction  f  /  *  * 

P3 


fiber  direction 


XFEM  model: 


active  phase 


Figure  2.73-1:  Composite  plate  model  for  PAC  validation  study. 
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Figure  2. 7. 3-2:  Comparison  of  experimental  and  XFEM  results:  fixity  over  fiber  orientation. 

programmed  stage  (ii)  it  is  clamped  at  its  center  to  suppress  rigid  body  motion.  In  the  case  of  the  single 
stage  model,  the  eigenstrain  is  assumed  to  be  £y  =  diag[Ul  / L,-0.5Ul  /  L,-0.5Ul  I L\,  which  is 

consistent  with  the  programming  displacements  and  incompressible  material.  Details  on  the  level  set 
function  that  described  the  fiber  layout,  the  XFEM  model,  and  the  material  properties  are  given  in  Maute 
et  al.  [30]. 

The  response  of  the  PAC  is  studied  for  7  different  fiber  orientations.  The  experimental  results  of 
KM-Refl  and  the  XFEM  predictions  are  compared  in  Figure  2. 7. 3-2.  The  fixity  of  the  composite  plate  is 
measured  by  the  elongation  of  the  plate  in  the  programmed  stage  (ii)  relative  to  prescribed  elongation 
U]  at  the  programming  stage  (i).  The  fixity  parameter  of  the  active  phase  is  calibrated  by  considering  the 

experimental  results  for  a  plate  made  of  the  active  phase  only;  see  the  dotted  line  in  Figure  2. 7. 3-2. 
Comparing  the  experimentally  measured  and  numerical  results  shows  that  the  two-stage  model  captures 
well  the  PAC  response  across  a  wide  range  of  fiber  orientations.  However,  the  single  stage  model  leads 

only  to  acceptable  results  for  fiber  angles  up  to  cp  —  20  ,  because  the  assumption  for  the  eigenstrain 
distribution  is  valid  only  for  small  differences  between  training  load  and  fiber  directions. 

Given  the  complexity  of  the  material  behavior  and  the  uncertainty  of  the  experimental  results[31], 
the  two-stage  model  and,  to  a  lesser  extent,  the  single  stage  model  showed  satisfactory  agreement  with 
the  experimental  measurements  for  the  example  considered  above.  This  motivated  the  application  of 
these  PAC  models  to  the  design  optimization  problems  presented  in  the  following  section.  However, 
additional  validation  problems  need  to  be  studied  in  the  future  to  increase  the  confidence  in  the  proposed 
PAC  model.  Future  studies  should  also  include  problems  that  exhibit  large  out-of-plane  deformations  and 
more  complex  inclusion  geometries. 


2.7.4  Optimization  of  PACs  by  LSM-XFEM  topology  optimization 

We  studied  the  proposed  LSM-XFEM  approach  for  PACs  with  several  design  problems,  using  the  single 
stage  and  two-stage  PAC  models.  Here  we  report  on  three  of  them.  In  the  first  two  examples,  the  target 
shape  promotes  a  plate  bending  type  response  of  increasing  complexity.  The  third  example  explores  the 
ability  to  generate  localized  deformations.  Measurements  for  3D  printed  realizations  of  the  optimized 
designs  were  compared  against  numerical  predictions  of  the  programmed  shape  at  stage  (ii).  The  design 
studies  illustrated  the  ability  of  the  proposed  optimization  approach  to  find  non-trivial  designs.  The 
comparison  against  experimental  measurements  motivates  future  research  on  improving  the  accuracy  of 
the  PAC  model. 


60 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


top  view: 


a;2  Vx 


Xi 

side  view: 
X3  U\ 


Ur 


-^target 


top  view: 


L* 


J-S 

S)  pj 
*33 

o>  2 

^  3 

i  O 

O-  © 

O  V 

3  % 

O  ^ 

P  c§* 

Ur 


xr 

side  view: 
Ui  x3  Ux 


passive  phase  p  =  1 


xr 


passive  phase  p  =  1  active  phase  p  =  2 

(a) 


-t 


ptarpet 


M 

•••••• 


p  targy 


Ux 


xr 


T  active  phase  p  =  2 

(b) 


Figure  2. 7. 4-1.  (a)  Setup  for  3D  PAC  layout  optimization  problem  --  plate  bending  configuration,  (b)  Setup 
for  3D  PAC  layout  optimization  problem  -  localized  deformation  configuration. 


For  all  examples,  we  study  the  formulation  of  the  optimization  problem:  we  minimize  the  square  of  the 
difference  between  the  displacements  in  the  programmed  stage  (ii)},  ui ,  and  a  target  deformation  ufxgQt 

integrated  over  the  surface  rtarget  where  the  deformations  are  monitored.  This  yields  the  following 
contribution  to  the  objective  function: 

dT •  i2-7-4-1' 

where  wt  denotes  the  weighting  factors  for  the  components  of  the  displacement  vector.  To  prevent  the 

formation  of  small  geometric  features  and  irregular  shapes  we  augment  the  objective  ^rtaiget  with  the 
following  regularization  term  that  penalizes  the  perimeter  of  the  phase  boundary: 

f"‘=> vfr  dT,  12.7.4-2) 

1,2 

where  y  denotes  the  perimeter  penalty  factor.  The  value  of  the  perimeter  penalty  factor  is  chosen  to 
smooth  the  surface  of  the  phase  boundary  without  significantly  affecting  the  shape  matching  objective 
jFtarget .  The  choice  of  y  is  problem  dependent  and  has  been  determined  iteratively.  While  for  the 

problems  studied  penalizing  the  perimeter  often  led  to  satisfactory  results. 

The  geometry  of  the  samples  and  boundary  conditions  in  the  programming  stage  (i)  are  depicted 
in  Figures  2.7.4-la  and  b  for  the  plate  bending  and  for  the  localized  deformation  problems,  respectively. 
Details  on  the  dimensions  and  material  parameters  are  given  in  Maute  et  al. [30].  The  difference  between 
the  target  and  programmed  shape  is  measured  at  the  upper  surface  for  the  first  three  examples  and  at 
the  upper  and  lower  surfaces  for  the  fourth  problem.  The  analytical  expressions  for  the  target 
displacements  are  given  in  Table  2. 7.4-1  and  are  defined  with  respect  to  a  coordinate  system  located  at 
the  center  of  the  design  domain. 

In  the  programming  stage  (i),  the  samples  are  stretched  by  6%  by  applying  prescribed 
displacements  on  front  and  rear  faces.  In  the  programmed  stage  (ii),  the  samples  are  supported  at  their 
center  to  prevent  rigid  body  motion.  In  all  examples,  the  design  domain  is  discretized  by  8-node,  tri-linear 
finite  elements.  Details  on  the  mesh  sizes  and  solver  strategies  are  given  in  Maute  et  al.  [30]. 
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Table  2. 7.4-1.  Target  shapes  for  PAC  layout  problems. 


Case 

u3 ,  mm 

Parabolic  bending 

-2.4 

fx0 

Uoj 

2 

Twisted  parabolic  bending 

-4.8 

M 

Uoj 

2  s  \ 

.  7TX1  1 

sm  — - 

Uoj 

Localized  deformation 

-0.1  sign(x3 )  max(0, 1  -  (xf  +  x\ )  /  25) 

The  level  set  function  is  discretized  by  the  same  mesh  used  for  the  PAC  analysis.  The  nodal  level 
set  values  are  defined  as  functions  of  the  optimization  variables,  using  a  linear  filter[30].  The  resulting 
parameter  optimization  problems  are  solved  by  the  Globally  Convergent  Method  of  Moving 
Asymptotes[32]  without  sub-cycling.  The  design  sensitivities  of  the  objective  and  constraints  are 
determined  by  the  adjoint  method.  The  designs  are  initialized  with  uniform  arrays  of  cuboids  of  either 
active  or  passive  material,  surrounded  by  material  of  the  other  phase. 

The  optimized  designs  were  3D  printed  by  exporting  the  triangulated  phase  boundaries  of  the 
XFEM  models  for  each  phase  separately,  using  the  STereo  Lithography  (STL)  file  format.  To  this  end,  we 
extracted  the  XFEM  mesh  used  for  integrating  the  weak  form  of  the  governing  equations.  Note  that  this 
triangulation  is  aligned  with  the  phase  boundaries  and  thus  each  STL  file  defines  the  volume  of  one  distinct 
phase.  The  STL  files  were  loaded  into  the  software  interface  of  an  Objet  Connex  260  printer  and  a  material 
type  was  assigned  to  each  volume  associated  with  an  STL  file.  The  printer  software  optimizes  the  material 
layout  along  the  phase  boundaries  for  improved  bonding  and  material  strength.  The  printed  PAC  sample 
were  subject  to  the  thermo-mechanical  programming  cycle  described  above  and  the  programmed  shape 
in  stage  (ii)  was  measured.  Note  that  the  experimental  results  reported  here  mainly  serve  as  a  proof  of 
concept  but  lack  the  accuracy  needed  for  rigorous  quantitative  comparisons  between  predictions  and 
experimental  measurements. 

In  Figures  2. 7. 4-2,  Figure  2. 7. 4-3,  and  Figure  2. 7.4-4  the  results  for  the  three  problems  are  shown. 
In  each  figure  the  topology  optimization  result  for  the  single  and  two-stage  models  are  shown  with  the 
active  phase  colored  in  dark  gray;  the  bright  gray  color  is  the  inactive  phase.  The  figures  also  show  the 
finite  element  prediction  of  the  two-stage  optimized  design  with  the  colors  representing  the  u3 

displacement  component.  The  printed  design  is  shown  in  the  second  row.  Note  for  visualization  purposes, 
we  print  the  sample  with  a  transparent  matrix  material  and  a  dark  SMP  material.  However,  the  mechanical 
test  is  performed  with  Tango  Black+  and  Shore95.  The  deformed  sample  after  activation  is  shown  for  the 
first  two  examples.  However,  as  our  experimental  instrumentation  lacked  the  ability  to  sufficiently  resolve 
sub-millimeter  deformations,  we  were  not  able  to  obtain  meaningful  measurements  for  the  third 
problem. 


Single  stage  design: 


T wo-stage  design : 


Two-stage  design: 


u.3  \rmn] 
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Sample  3D  printed  with  transparent  passive  phase: 
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Figure  2. 7.4-2:  Optimization  results  for  paraboli 


Programmed  sample  at  stage  (ii): 


bending  target  shape. 


Single  stage  design:  Two-stage  design: 


Two-stage  design: 


Figure  2. 7.4-3:  Optimization  results  for  twisted  parabolic  bending  target  shape. 
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The  numerical  results  demonstrated  that  the  proposed  optimization  approach  provides  a  crisp 
and  well  resolved  description  of  the  phase  boundaries.  Surface  meshes  of  the  phase  boundaries  could  be 
easily  extracted  from  the  XFEM  model  and  directly  exported  into  the  software  interface  of  3D  printers 
without  the  need  for  any  form  of  post-processing. 

For  all  design  problems  studied  with  a  bending  type  target  deformation,  the  target  shapes  could 
be  matched  with  acceptable  accuracy.  Both  PAC  models  resulted  in  conceptually  similar  designs  with 
sophisticated  spatial  arrangements  of  the  active  phase.  The  results  for  the  example  with  a  localized  target 
deformation  suggest  that  creating  localized  deformations  is  a  challenging  design  problem  even  if  the 
amplitude  of  the  target  design  is  small.  Nevertheless,  the  examples  studied  to  date  demonstrate  the 
ability  of  the  proposed  optimization  method  to  find  non-trivial  and  non-intuitive  designs.  Analyzing  the 
designs  optimized  with  the  single  stage  PAC  model  by  the  more  accurate  two-stage  PAC  model  showed 
that  the  shape  mismatch  is  significantly  larger  than  the  one  of  the  designs  optimized  with  the  two-stage 
model,  see[30].  This  observation  suggests  that  the  accuracy  of  the  PAC  model  has  a  noticeable  influence 
on  the  optimization  results  and  therefore  an  accurate  PAC  model  should  be  used  for  topology 
optimization. 

The  comparison  of  the  XFEM  predictions  against  experimental  measurements  for  the  optimized 
designs  of  the  examples  with  a  bending  type  target  deformation  suggests  that  the  proposed  two-stage 
PAC  model  may  be  appropriate  for  design  purposes  but  insufficient  to  accurately  predict  the  PAC  response 
in  the  presence  of  large  deformations.  During  the  third  year  of  this  project,  more  accurate  but  still 
computationally  tractable  PAC  models  will  be  researched  and  integrated  into  the  proposed  optimization 
framework.  Furthermore,  the  perimeter  constraint  was  found  to  be  insufficient  to  explicitly  control  the 
local  feature  size  of  the  active  phase.  Thus,  methods  will  be  explored  to  constrain  the  minimum  feature 
size  and  to  account  for  the  resolution  limits  of  3D  printers. 

2.7.5  Finite  Strain  PAC  Model  for  Design  Optimization 

The  small  strain  assumption  on  which  the  PAC  models  described  in  Section  2.7.3  are  based  is  not  valid  for 
the  level  of  training  loads  used  for  practical  applications  and  does  not  capture  accurately  the  deformations 
achieved  upon  activation.  Therefore,  we  expanded  our  simplified  PAC  model  onto  finite  strains  and 
studied  the  influence  of  the  finite  strain  model  on  the  optimum  material  layout.  This  study  is  currently 
completed  and  a  paper  on  the  results  of  these  studies  is  currently  in  preparation. 

Finite  strain  PAC  model:  While  there  are  several  constitutive  models  available  to  accurately  describe  the 
behavior  of  amorphous  shape  memory  polymers  (SMPs),  these  models  are  too  complex  and  too 
computationally  costly  for  design  optimization.  Previously,  we  have  created  an  analysis  approach  more 
amenable  for  design  in  terms  of  computational  speed  and  ease  of  implementation  but  with  acceptable 
fidelity.  This  approach  accounts  for  stress  relaxation  and  the  resulting  residual  stresses  via  eigenstrains 
that  can  vary  with  position.  The  eigenstrains  are  predicted  by  resolving  the  strain  field  in  the  PAC  at  the 
training  stage  due  to  the  programming  loads  at  temperatures  above  the  glass  transition  temperature. 
However,  our  previous  PAC  model  assumes  infinitesimal  strains  and  a  linear  material  behavior.  To 
overcome  this  limitation,  we  developed  a  nonlinear  PAC  model  which  accounts  for  finite  strains,  but  is 
still  considerably  simpler  than  constitutive  models  that  describe  the  phase  transition  and  the  associated 
changes  in  material  properties  during  the  shape  memory  cycle. 

We  consider  the  response  of  the  PAC  at  the  training  stage  subject  to  programming  loads,  referred 
to  stage  (i),  and  after  the  programming  loads  are  removed,  referred  to  stage  (ii).  The  effects  of  the  phase 
transition  of  the  active  phase  are  modeled  via  eigenstrains,  which  are  computed  in  stage  (i)  and  applied 
in  stage  (ii),  and  a  change  in  Young's  modulus.  The  residuals  of  the  static  equilibrium  equations  in  both 
stages  are  written  in  weak  form  using  a  total  Lagrangian  formulation  as  follows: 
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(2.7.5-1) 
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(2. 7. 5-3) 


where  (~)  denotes  quantities  in  stage  (i)  and  (~)  in  stage  (ii).  Passive  and  active  phases  are  denoted 
with  superscripts  ”1"  and  "2".  The  displacement  vector,  the  deformation  gradient  tensor,  and  the  1st 
Piola-Kirchhoff  stress  tensor  are  denoted  by  uf  ,  F and  /^^respectively.  External  forces,  T: p ,  act  along 

the  boundary  rp0  T  in  the  reference  configuration,  and  prescribed  displacements,  UF ,  are  applied  at  the 

boundary  .  Note  the  boundary  conditions  typically  differ  between  stages  (i)  and  (ii).  For  the  sake  of 
simplicity,  the  material  behavior  at  the  two  stages  is  described  by  a  Saint-Venant  Kirchhoff  model: 


S’ = 


(2.7.5-41 


where  SF  and  EF  denote  the  2nd  Piola-Kirchhoff  stress  tensor  and  the  Green-Lagrange  strain  tensor, 

respectively.  Other  hyperelastic  material  models  can  be  easily  accommodated.  The  Lame  constants  are 
A  and  /u  .  Green-Lagrange  strain  tensors  is  computed  for  follows: 


(2. 7. 5-5) 


where  FF'P  is  the  elastic  contribution  of  the  deformation  gradient  tensor.  In  the  training  stage,  we  assume 
that  there  are  no  inelastic  contributions  to  the  total  deformation  gradient  tensor,  Fp ,  such  that: 

Fr-F-^.  12.7.5-6, 

where  X .  denotes  the  location  of  a  material  point  in  the  reference  configuration.  In  the  actuation  stage 

(ii),  we  model  the  shape  memory  effect  via  a  multiplicative  decomposition  of  the  deformation  gradient 
tensor: 

F’ =F”fy’  with  Fp=SiJ+^-,  (2. 7. 5-7) 

where  the  inelastic  contribution,  FF'P ,  represents  the  eigenstrains  due  to  deformations  at  the 
programming  stage  and  is  modeled  as  follows: 

FF’p=fpFp,  (2.7 .5-8) 
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where  /  is  the  strain  fixity  parameter.  Note  that  this  approach  is  similar  to  the  proposed  in  [34]  but  does 
not  require  modeling  the  phase  transition.  Instead,  the  fixity  parameter  in  (8)  represents  the  complete 
stress  relaxation  process  and  needs  to  be  calibrated  for  particular  thermos-mechanical  cycling  processes. 
The  Young's  modulus  of  the  active  phase  increases  as  the  material  transitions  from  a  rubbery  state  in 
stage  (i)  to  a  glassy  state  in  stage  (ii).  The  elastic  tensor  of  the  passive  phase  is  assumed  constant.  For  the 

passive  phase  the  fixity  is  assumed  negligible,  i.e.  /’  =  0.0  . 


XFEM  discretization  of  PAC  model:  The  above  finite  strain  PAC  model  was  discretized  by  a  Heaviside 
enriched  formulation  of  the  XFEM.  In  previous  studies,  we  have  shown  that  this  enrichment  strategy 
provides  great  flexibility  and  accuracy  in  modeling  heterogeneous  materials  [27,  28].  However,  as  the 
Heaviside  enrichment  does  not  intrinsically  guarantee  the  continuity  of  the  displacement  fields  across 
material  interfaces,  additional  continuity  constraints  need  to  be  introduced.  In  this  work,  we  use  a  Nitsche 
formulation  [35].  Furthermore,  the  XFEM  problem  may  become  ill-conditioned  when  elements  are 
interested  such  that  small  volumes  of  one  phase  are  generated.  These  configurations  do  not  only  cause 
issues  for  solving  the  set  of  nonlinear  equations,  e.g.  by  Newton's  method,  but  also  lead  to  stress 
oscillations  along  the  interface.  The  latter  may  severely  affect  the  stability  of  the  XFEM  formulation  of 
finite  strain  problems.  To  overcome  these  stability  issues,  we  adopted  and  extended  the  so-called  ghost 
penalization  approach  of  Burman  and  Hansbo  [36,  37]  for  finite  strain  problems. 

The  weak  form  of  the  governing  equations  augmented  by  Nitsche's  formulation  of  the  interface 
continuity  condition  and  ghost  stabilization  can  be  written  as  follows: 


+  Z  X  Ygp\f 


P=1  F»sFP 


f 

s 

'  du?' 

n 

j  duf  "1 

n. 

Jf/ 

w 

j 

K 

dF 


(2. 7. 5-9) 


■  Jr„  [<*“. ] {Pfa )  dr + Jr„  [8P'«, ] {«, }  dT  +  2  fK k  [Su,  ] [gu,  ]  dr, 


e&Er. 


where  the  second  row  is  the  contribution  of  the  ghost  penalty  formulation  and  the  terms  in  the  third  row 
are  due  to  the  Nitsche  interface  condition.  The  symbol  [•] denotes  the  jump  across  an  interface;  the 

symbol  {•}  represents  an  averaged  value. 

The  ghost  penalization  measures  the  jump  of  the  displacement  gradients  across  an  element 
interfaces,  projected  onto  the  interface  normal.  The  set  Fcput  includes  all  interfaces  between  elements 
where  at  least  one  of  the  elements  is  intersected  and  the  interface  contains  the  phase  p,  i.e. 
Fp  ^0.  Note  that  these  integrals  are  independent  of  the  intersection  configuration  but  are 

evaluated  across  the  entire  element  fac eFp ,  such  that  the  degrees  of  freedom  interpolating  the 
displacement  field  up  maintain  a  non-vanishing  zone  of  influence.  The  penalty  parameter  yGP  controls 

the  strength  of  the  ghost  penalty  term.  The  larger  yGP,  the  smoother  the  spatial  gradients  of  the 

displacement  field  and,  thus,  the  stress  field  along  the  material  interfaces,  increasing  the  stability  of  the 
XFEM  formulation  at  the  cost  of  accuracy. 


Design  optimization  study  with  nonlinear  of  PAC  model:  To  investigate  the  influence  of  accounting  for 
finite  strains  on  the  optimal  PAC  layout,  we  revisited  a  design  study  we  performed  earlier  using  our  linear 
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PAC  model  [30].  The  design  problem  is  to  find  the  arrangement  of  the  SMP  phase  within  a  rubbery  passive 
matrix  such  that  a  plate  assumes  a  target  shape  after  activation.  To  this  end,  we  minimize  the  square  of 
the  difference  between  the  displacements  in  the  programmed  stage  (ii)},  ui ,  and  a  target  deformation 

i/*arget  integrated  over  the  surface  rtaiget  where  the  deformations  are  monitored.  This  yields  the  following 
contribution  to  the  objective  function: 

=  Jr„  ",  («,  -  f  dr,  (2.7.5-10) 


where  wt  denotes  the  weighting  factors  for  the  components  of  the  displacement  vector.  To  prevent  the 

formation  of  small  geometric  features  and  irregular  shapes  we  augment  the  objective  ^?target  with  the 
following  regularization  term  that  penalizes  the  perimeter  of  the  phase  boundary: 

•^reg=7„Jr  dr,  (2.7.5-11) 

1,2 


where  y  denotes  the  perimeter  penalty  factor.  The  value  of  the  perimeter  penalty  factor  is  chosen  to 
smooth  the  surface  of  the  phase  boundary  without  significantly  affecting  the  shape  matching  objective 

^-target  c|10jce  0f 

yper  is  problem  dependent  and  has  been  determined  iteratively.  While  for  the 

problems  studied  penalizing  the  perimeter  often  led  to  satisfactory  results. 

The  geometry  of  the  plate  and  the  boundary  conditions  in  the  programming  stage  (i)  are  depicted 
in  Figures  2. 7. 5-1.  The  target  shape  represents  a  parabola  with: 


^  target 


2.4 


Xl 

20 


w1  =w2  =0 


w3=  1 


(2.7.5-12) 


In  the  programming  stage  (i),  the  samples  are  stretched  by  6%  by  applying  prescribed  displacements  on 
front  and  rear  faces.  In  the  programmed  stage  (ii),  the  samples  are  supported  at  their  center  to  prevent 
rigid  body  motion.  In  all  examples,  the  design  domain  is  discretized  by  8-node,  tri-linear  finite  elements. 
The  forward  problem  is  solved  by  a  nonlinear  block  Gauss-Seidel  algorithm,  thus  taking  advantage  that 
that  the  stage  (i)  problem  is  independent  of  the  stage  (ii)  problem.  We  observed  that  this  strategy  provides 
significant  increase  in  robustness  and  reduction  in  computational  costs  over  solving  simultaneously  the 
stage  (i)  and  (ii)  problems. 


top  view: 


Ui 
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Figure  2. 7. 5-1:  Geometry  and  boundary  conditions  of  PAC  design  optimization  problem. 
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Figure  2. 7. 5-2:  Optimization  results  using  finite  strain  PAC  models. 

The  level  set  function  is  discretized  by  the  same  mesh  used  for  the  PAC  analysis.  The  nodal  level  set  values 
are  defined  as  functions  of  the  optimization  variables,  using  a  linear  filter  [30].  The  resulting  parameter 
optimization  problems  are  solved  by  the  Globally  Convergent  Method  of  Moving  Asymptotes[32]  without 
sub-cycling.  The  design  sensitivities  of  the  objective  and  constraints  are  determined  by  the  adjoint 
method.  Here  we  solve  the  adjoint  equations  of  the  stage  (i)  and  (ii)  problems  simultaneously,  as  these 
problems  are  linear.  The  designs  are  initialized  with  uniform  arrays  of  cuboids  of  either  active  or  passive 
material,  surrounded  by  material  of  the  other  phase. 

Figure  2. 7. 5-2  shows  the  layout  of  the  SMP  phase  in  the  optimized  design  obtained  with  the 
nonlinear  PAC  model.  Comparing  this  layout  with  the  one  obtained  previously  for  a  small  strain  PAC  model 
shows  that  the  designs  have  a  similar  geometry  but  differ  locally.  Note  the  additional  ribs  obtained  when 
using  the  nonlinear  model.  The  design  optimized  with  the  linear  PAC  model  is  inferior  to  the  one  optimized 
with  the  nonlinear  PAC  model.  In  future  studies,  additional  optimization  problems  using  the  nonlinear 
PAC  model  will  be  performed  and  the  potential  of  utilizing  large  deformation  and  instability  phenomena, 
such  as  buckling,  will  be  investigated. 

2.7.6  Level-Set  Topology  Optimization  with  Contact,  Cohesion,  and  Debonding 

An  important  design  consideration  for  engineered  multiphase  materials  is  the  mechanical  behavior  along 
the  material  interfaces,  in  particular  contact,  cohesions,  and  debonding.  To  account  for  these  phenomena 
in  the  design  optimization  process,  we  expanded  our  unique  LSM-XFEM  design  optimization  framework 
to  contact  and  cohesive  interface  problems.  In  contrast  to  density  based  topology  optimization 
approaches,  where  such  phenomena  cannot  be  resolved,  the  crispness  of  the  interface  description  in  both 
the  geometry  and  analysis  models  allows  us  to  solve  and  study  such  problems.  The  following  study  is 
published  in:  "R.  Behrou,  M.  Lawry,  and  K.  Maute.  Level  set  topology  optimization  of  structural  problems 
with  interface  cohesion,  Int.  J.  Num.  Meth.  Engrg,  2016,  submitted". 

Topology  Optimization  considering  cohesive  interface  behavior:  Building  upon  the  LSM-XFEM 
framework  developed  in  previous  years  as  part  of  this  project,  we  created  a  finite  element  topology 
optimization  approach  for  the  design  of  two-phase  structural  systems  considering  contact  and  cohesion 
phenomena  along  material  interfaces.  The  geometry  of  the  material  interface  is  described  by  an  explicit 
level  set  method  and  the  structural  response  is  predicted  by  the  extended  finite  element  method.  In  this 
work  the  interface  condition  is  described  by  a  bilinear  cohesive  zone  model  based  on  the  traction- 
separation  constitutive  relation.  The  non-penetration  condition  in  the  presence  of  compressive  interface 
forces  is  enforced  by  a  stabilized  Lagrange  multiplier  method.  The  mechanical  model  assumes  a  linear 
elastic  isotropic  material,  infinitesimal  strain  theory,  and  a  quasi-static  response.  The  optimization 
problem  is  solved  by  a  nonlinear  programming  method  and  the  design  sensitivities  are  computed  by  the 
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adjoint  method.  For  details  on  the  mechanical  model  and  the  XFEM  formulation  of  the  cohesive  interface 
model,  the  reader  is  referred  to  [38]. 

The  performance  of  the  method  outlined  above  was  evaluated  studied  with  the  optimization 
problem  depicted  in  Figure  2. 7. 6-1,  considering  both  a  2D  and  a  3D  configuration.  A  structural  anchor 
(represented  by  phase  A)  is  embedded  in  a  host  material  (designated  as  phase  B)  with  frictionless  and 
cohesive  interface  conditions,  imposed  at  the  boundary  between  both  phases.  We  wish  to  determine  the 
optimal  geometry  such  that  the  holding  force  of  the  anchor  is  maximized. 


Figure  2. 7. 6-1:  Design  problems  with  cohesive  contact  interfaces. 
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Figure  2. 7. 6-2:  Evolution  of  2D  design  assuming  frictionless  contact  along  the  material  interfaces. 
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Figure  2. 7. 6-3:  Influence  of  interface  conditions  on  the  optimized  design  (2D  case). 
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Figure  2. 7. 6-4:  Evolution  of  3D  design  assuming  frictionless  contact  along  the  material  interfaces. 
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Figure  2. 7. 6-5:  Influence  of  interface  conditions  on  the  optimized  design  (3D  case). 

Figure  2. 7. 6-2  shows  the  evolution  of  the  design  for  the  2D  case,  assuming  no-cohesion  but 
accounting  for  frictionless  contact.  The  holding  force  increases  by  forming  noticeable  barb-like  structures. 
In  Figure  2. 7. 6-3  we  compare  the  influence  of  assuming  cohesive  forces  in  normal  and  tangential 
directions.  We  vary  the  resistance  in  horizontal  direction  along  the  horizontal  axis  and  the  resistance  in 
normal  direction  along  the  vertical  axis.  The  result  in  the  lower  left  corner  is  the  one  for  frictionless  contact 
and  the  result  in  the  upper  right  corner  for  perfect  bonding.  This  design  study  illustrates  the  significant 
impact  of  the  material  interface  properties  on  the  optimal  design. 

In  Figure  2. 7. 6-4  and  Figure  2. 7. 6-5,  the  results  for  the  3D  configuration  are  shown.  Similar  to  the 
2D  design  with  frictionless  contact,  also  the  optimized  3D  design  features  barbs  along  the  interface.  Again, 
the  interface  conditions  have  a  significant  influence  on  the  optimized  design.  More  details  can  be  found 
in  [38]. 

In  summary,  the  results  obtained  from  topology  optimization  reveal  distinct  design  characteristics 
for  the  various  interface  phenomena  considered.  In  addition,  three-dimensional  examples  demonstrate 
optimal  geometries  that  cannot  be  fully  captured  by  reduced  dimensionality.  The  optimization  framework 
presented  is  limited  to  two-phase  structural  systems  where  the  material  interface  is  coincident  in  the 
undeformed  configuration,  and  to  structural  responses  that  remain  valid  considering  small  strain 
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kinematics.  In  an  ongoing  study,  the  above  results  have  been  expanded  onto  accounting  for  finite  strains 
[39]. 
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